How to determine if a list of polygon points are in clockwise order?

Having a list of points, how do I find if they are in clockwise order?

For example:

point[0] = (5,0)
point[1] = (6,4)
point[2] = (4,5)
point[3] = (1,5)
point[4] = (1,0)

would say that it is anti-clockwise (or counter-clockwise, for some people).

Some of the suggested methods will fail in the case of a non-convex polygon, such as a crescent. Here's a simple one that will work with non-convex polygons (it'll even work with a self-intersecting polygon like a figure-eight, telling you whether it's mostly clockwise).

Sum over the edges, (x₂ − x₁)(y₂ + y₁). If the result is positive the curve is clockwise, if it's negative the curve is counter-clockwise. (The result is twice the enclosed area, with a +/- convention.)

point[0] = (5,0)   edge[0]: (6-5)(4+0) =   4
point[1] = (6,4)   edge[1]: (4-6)(5+4) = -18
point[2] = (4,5)   edge[2]: (1-4)(5+5) = -30
point[3] = (1,5)   edge[3]: (1-1)(0+5) =   0
point[4] = (1,0)   edge[4]: (5-1)(0+0) =   0
                                         ---
                                         -44  counter-clockwise

Find the vertex with smallest y (and largest x if there are ties). Let the vertex be A and the previous vertex in the list be B and the next vertex in the list be C. Now compute the sign of the cross product of AB and AC.

---

References:

• How do I find the orientation of a simple polygon? in Frequently Asked Questions: comp.graphics.algorithms.

• Curve orientation at Wikipedia.

I'm going to throw out another solution because it's straightforward and not mathematically intensive - it just uses basic algebra. Calculate the signed area of the polygon. If it's negative the points are in clockwise order, if it's positive they are counterclockwise. (This is very similar to Beta's solution.)

Calculate the signed area: A = 1/2 * (x₁*y₂ - x₂*y₁ + x₂*y₃ - x₃*y₂ + ... + x_n*y₁ - x₁*y_n)

Or in pseudo-code:

signedArea = 0
for each point in points:
    x1 = point[0]
    y1 = point[1]
    if point is last point
        x2 = firstPoint[0]
        y2 = firstPoint[1]
    else
        x2 = nextPoint[0]
        y2 = nextPoint[1]
    end if

    signedArea += (x1 * y2 - x2 * y1)
end for
return signedArea / 2

Note that if you are only checking the ordering, you don't need to bother dividing by 2.

Sources: http://mathworld.wolfram.com/PolygonArea.html

The cross product measures the degree of perpendicular-ness of two vectors. Imagine that each edge of your polygon is a vector in the x-y plane of a three-dimensional (3-D) xyz space. Then the cross product of two successive edges is a vector in the z-direction, (positive z-direction if the second segment is clockwise, minus z-direction if it's counter-clockwise). The magnitude of this vector is proportional to the sine of the angle between the two original edges, so it reaches a maximum when they are perpendicular, and tapers off to disappear when the edges are collinear (parallel).

So, for each vertex (point) of the polygon, calculate the cross-product magnitude of the two adjoining edges:

Using your data:
point[0] = (5, 0)
point[1] = (6, 4)
point[2] = (4, 5)
point[3] = (1, 5)
point[4] = (1, 0)

So Label the edges consecutively as
edgeA is the segment from point0 to point1 and
edgeB between point1 to point2
...
edgeE is between point4 and point0.

Then Vertex A (point0) is between
edgeE [From point4 to point0]
edgeA [From point0 to `point1'

These two edges are themselves vectors, whose x and y coordinates can be determined by subtracting the coordinates of their start and end points:

edgeE = point0 - point4 = (1, 0) - (5, 0) = (-4, 0) and
edgeA = point1 - point0 = (6, 4) - (1, 0) = (5, 4) and

And the cross product of these two adjoining edges is calculated using the determinant of the following matrix, which is constructed by putting the coordinates of the two vectors below the symbols representing the three coordinate axis (i, j, & k). The third (zero)-valued coordinate is there because the cross product concept is a 3-D construct, and so we extend these 2-D vectors into 3-D in order to apply the cross-product:

 i    j    k 
-4    0    0
 1    4    0    

Given that all cross-products produce a vector perpendicular to the plane of two vectors being multiplied, the determinant of the matrix above only has a k, (or z-axis) component.
The formula for calculating the magnitude of the k or z-axis component is
a1*b2 - a2*b1 = -4* 4 - 0* 1 = -16

The magnitude of this value (-16), is a measure of the sine of the angle between the 2 original vectors, multiplied by the product of the magnitudes of the 2 vectors.
Actually, another formula for its value is
A X B (Cross Product) = |A| * |B| * sin(AB).

So, to get back to just a measure of the angle you need to divide this value, (-16), by the product of the magnitudes of the two vectors.

|A| * |B| = 4 * Sqrt(17) = 16.4924...

So the measure of sin(AB) = -16 / 16.4924 = -.97014...

This is a measure of whether the next segment after the vertex has bent to the left or right, and by how much. There is no need to take arc-sine. All we will care about is its magnitude, and of course its sign (positive or negative)!

Do this for each of the other 4 points around the closed path, and add up the values from this calculation at each vertex..

If final sum is positive, you went clockwise, negative, counterclockwise.

Here is a simple C# implementation of the algorithm based on @Beta's answer.

Let's assume that we have a Vector type having X and Y properties of type double.

public bool IsClockwise(IList<Vector> vertices)
{
    double sum = 0.0;
    for (int i = 0; i < vertices.Count; i++) {
        Vector v1 = vertices[i];
        Vector v2 = vertices[(i + 1) % vertices.Count];
        sum += (v2.X - v1.X) * (v2.Y + v1.Y);
    }
    return sum > 0.0;
}

% is the modulo or remainder operator performing the modulo operation which (according to Wikipedia) finds the remainder after division of one number by another.

---

Optimized version according to @MichelRouzic's comment:

double sum = 0.0;
Vector v1 = vertices[vertices.Count - 1]; // or vertices[^1] with
                                          // C# 8.0+ and .NET Core
for (int i = 0; i < vertices.Count; i++) {
    Vector v2 = vertices[i];
    sum += (v2.X - v1.X) * (v2.Y + v1.Y);
    v1 = v2;
}
return sum > 0.0;

This saves not only the modulo operation % but also an array indexing.

An implementation of Sean's answer in JavaScript:

function calcArea(poly) {
    if(!poly || poly.length < 3) return null;
    let end = poly.length - 1;
    let sum = poly[end][0]*poly[0][1] - poly[0][0]*poly[end][1];
    for(let i=0; i<end; ++i) {
        const n=i+1;
        sum += poly[i][0]*poly[n][1] - poly[n][0]*poly[i][1];
    }
    return sum;
}

function isClockwise(poly) {
    return calcArea(poly) > 0;
}

let poly = [[352,168],[305,208],[312,256],[366,287],[434,248],[416,186]];

console.log(isClockwise(poly));

let poly2 = [[618,186],[650,170],[701,179],[716,207],[708,247],[666,259],[637,246],[615,219]];

console.log(isClockwise(poly2));

Pretty sure this is right. It seems to be working :-)

Those polygons look like this, if you're wondering:

Start at one of the vertices, and compute the angle subtended by each side.

The first and the last will be zero (so skip those); for the rest, the sine of the angle will be given by the cross product of the normalizations to unit length of (point[n]-point[0]) and (point[n-1]-point[0]).

If the sum of the values is positive, then your polygon is drawn in the anti-clockwise sense.

For what it is worth, I used this mixin to calculate the winding order for Google Maps API v3 apps.

The code leverages the side effect of polygon areas: a clockwise winding order of vertexes yields a positive area, while a counter-clockwise winding order of the same vertexes produces the same area as a negative value. The code also uses a sort of private API in the Google Maps geometry library. I felt comfortable using it - use at your own risk.

Sample usage:

var myPolygon = new google.maps.Polygon({/*options*/});
var isCW = myPolygon.isPathClockwise();

Full example with unit tests @ http://jsfiddle.net/stevejansen/bq2ec/

/** Mixin to extend the behavior of the Google Maps JS API Polygon type
 *  to determine if a polygon path has clockwise of counter-clockwise winding order.
 *  
 *  Tested against v3.14 of the GMaps API.
 *
 *  @author  stevejansen_github@icloud.com
 *
 *  @license http://opensource.org/licenses/MIT
 *
 *  @version 1.0
 *
 *  @mixin
 *  
 *  @param {(number|Array|google.maps.MVCArray)} [path] - an optional polygon path; defaults to the first path of the polygon
 *  @returns {boolean} true if the path is clockwise; false if the path is counter-clockwise
 */
(function() {
  var category = 'google.maps.Polygon.isPathClockwise';
     // check that the GMaps API was already loaded
  if (null == google || null == google.maps || null == google.maps.Polygon) {
    console.error(category, 'Google Maps API not found');
    return;
  }
  if (typeof(google.maps.geometry.spherical.computeArea) !== 'function') {
    console.error(category, 'Google Maps geometry library not found');
    return;
  }

  if (typeof(google.maps.geometry.spherical.computeSignedArea) !== 'function') {
    console.error(category, 'Google Maps geometry library private function computeSignedArea() is missing; this may break this mixin');
  }

  function isPathClockwise(path) {
    var self = this,
        isCounterClockwise;

    if (null === path)
      throw new Error('Path is optional, but cannot be null');
      
    // default to the first path
    if (arguments.length === 0)
        path = self.getPath();
 
    // support for passing an index number to a path
    if (typeof(path) === 'number')
        path = self.getPaths().getAt(path);

    if (!path instanceof Array && !path instanceof google.maps.MVCArray)
      throw new Error('Path must be an Array or MVCArray');
      
    // negative polygon areas have counter-clockwise paths
    isCounterClockwise = (google.maps.geometry.spherical.computeSignedArea(path) < 0);

    return (!isCounterClockwise);
  }

  if (typeof(google.maps.Polygon.prototype.isPathClockwise) !== 'function') {
    google.maps.Polygon.prototype.isPathClockwise = isPathClockwise;
  }
})();

This is the implemented function for OpenLayers 2. The condition for having a clockwise polygon is area < 0, it confirmed by this reference.

function IsClockwise(feature)
{
    if(feature.geometry == null)
        return -1;

    var vertices = feature.geometry.getVertices();
    var area = 0;

    for (var i = 0; i < (vertices.length); i++) {
	    j = (i + 1) % vertices.length;

	    area += vertices[i].x * vertices[j].y;
	    area -= vertices[j].x * vertices[i].y;
	    // console.log(area);
    }

    return (area < 0);
}

C# code to implement lhf's answer:

// https://en.wikipedia.org/wiki/Curve_orientation#Orientation_of_a_simple_polygon
public static WindingOrder DetermineWindingOrder(IList<Vector2> vertices)
{
	int nVerts = vertices.Count;
	// If vertices duplicates first as last to represent closed polygon,
	// skip last.
	Vector2 lastV = vertices[nVerts - 1];
	if (lastV.Equals(vertices[0]))
		nVerts -= 1;
	int iMinVertex = FindCornerVertex(vertices);
	// Orientation matrix:
	//     [ 1  xa  ya ]
	// O = | 1  xb  yb |
	//     [ 1  xc  yc ]
	Vector2 a = vertices[WrapAt(iMinVertex - 1, nVerts)];
	Vector2 b = vertices[iMinVertex];
	Vector2 c = vertices[WrapAt(iMinVertex + 1, nVerts)];
	// determinant(O) = (xb*yc + xa*yb + ya*xc) - (ya*xb + yb*xc + xa*yc)
	double detOrient = (b.X * c.Y + a.X * b.Y + a.Y * c.X) - (a.Y * b.X + b.Y * c.X + a.X * c.Y);

	// TBD: check for "==0", in which case is not defined?
	// Can that happen?  Do we need to check other vertices / eliminate duplicate vertices?
	WindingOrder result = detOrient > 0
			? WindingOrder.Clockwise
			: WindingOrder.CounterClockwise;
	return result;
}

public enum WindingOrder
{
	Clockwise,
	CounterClockwise
}

// Find vertex along one edge of bounding box.
// In this case, we find smallest y; in case of tie also smallest x.
private static int FindCornerVertex(IList<Vector2> vertices)
{
	int iMinVertex = -1;
	float minY = float.MaxValue;
	float minXAtMinY = float.MaxValue;
	for (int i = 0; i < vertices.Count; i++)
	{
		Vector2 vert = vertices[i];
		float y = vert.Y;
		if (y > minY)
			continue;
		if (y == minY)
			if (vert.X >= minXAtMinY)
				continue;

		// Minimum so far.
		iMinVertex = i;
		minY = y;
		minXAtMinY = vert.X;
	}

	return iMinVertex;
}

// Return value in (0..n-1).
// Works for i in (-n..+infinity).
// If need to allow more negative values, need more complex formula.
private static int WrapAt(int i, int n)
{
	// "+n": Moves (-n..) up to (0..).
	return (i + n) % n;
}

If you use Matlab, the function ispolycw returns true if the polygon vertices are in clockwise order.

As also explained in this Wikipedia article Curve orientation, given 3 points p, q and r on the plane (i.e. with x and y coordinates), you can calculate the sign of the following determinant

If the determinant is negative (i.e. Orient(p, q, r) < 0), then the polygon is oriented clockwise (CW). If the determinant is positive (i.e. Orient(p, q, r) > 0), the polygon is oriented counterclockwise (CCW). The determinant is zero (i.e. Orient(p, q, r) == 0) if points p, q and r are collinear.

In the formula above, we prepend the ones in front of the coordinates of p, q and r because we are using homogeneous coordinates.

This is my solution using the explanations in the other answers:

def segments(poly):
    """A sequence of (x,y) numeric coordinates pairs """
    return zip(poly, poly[1:] + [poly[0]])

def check_clockwise(poly):
    clockwise = False
    if (sum(x0*y1 - x1*y0 for ((x0, y0), (x1, y1)) in segments(poly))) < 0:
        clockwise = not clockwise
    return clockwise

poly = [(2,2),(6,2),(6,6),(2,6)]
check_clockwise(poly)
False

poly = [(2, 6), (6, 6), (6, 2), (2, 2)]
check_clockwise(poly)
True

I think in order for some points to be given clockwise all edges need to be positive not only the sum of edges. If one edge is negative than at least 3 points are given counter-clockwise.

My C# / LINQ solution is based on the cross product advice of @charlesbretana is below. You can specify a reference normal for the winding. It should work as long as the curve is mostly in the plane defined by the up vector.

using System.Collections.Generic;
using System.Linq;
using System.Numerics;

namespace SolidworksAddinFramework.Geometry
{
    public static class PlanePolygon
    {
        /// <summary>
        /// Assumes that polygon is closed, ie first and last points are the same
        /// </summary>
       public static bool Orientation
           (this IEnumerable<Vector3> polygon, Vector3 up)
        {
            var sum = polygon
                .Buffer(2, 1) // from Interactive Extensions Nuget Pkg
                .Where(b => b.Count == 2)
                .Aggregate
                  ( Vector3.Zero
                  , (p, b) => p + Vector3.Cross(b[0], b[1])
                                  /b[0].Length()/b[1].Length());

            return Vector3.Dot(up, sum) > 0;

        } 
        
    }
}

with a unit test

namespace SolidworksAddinFramework.Spec.Geometry
{
    public class PlanePolygonSpec
    {
        [Fact]
        public void OrientationShouldWork()
        {

            var points = Sequences.LinSpace(0, Math.PI*2, 100)
                .Select(t => new Vector3((float) Math.Cos(t), (float) Math.Sin(t), 0))
                .ToList();

            points.Orientation(Vector3.UnitZ).Should().BeTrue();
            points.Reverse();
            points.Orientation(Vector3.UnitZ).Should().BeFalse();



        } 
    }
}

After testing several unreliable implementations, the algorithm that provided satisfactory results regarding the CW/CCW orientation out of the box was the one, posted by OP in this thread (shoelace_formula_3).

As always, a positive number represents a CW orientation, whereas a negative number CCW.

A much computationally simpler method, if you already know a point inside the polygon:

1. Choose any line segment from the original polygon, points and their coordinates in that order.

2. Add a known "inside" point, and form a triangle.

3. Calculate CW or CCW as suggested here with those three points.

Here's swift 3.0 solution based on answers above:

    for (i, point) in allPoints.enumerated() {
        let nextPoint = i == allPoints.count - 1 ? allPoints[0] : allPoints[i+1]
        signedArea += (point.x * nextPoint.y - nextPoint.x * point.y)
    }
    
    let clockwise  = signedArea < 0

Another solution for this;

const isClockwise = (vertices=[]) => {
    const len = vertices.length;
    const sum = vertices.map(({x, y}, index) => {
        let nextIndex = index + 1;
        if (nextIndex === len) nextIndex = 0;
        
        return {
            x1: x,
            x2: vertices[nextIndex].x,
            y1: x,
            y2: vertices[nextIndex].x
        }
    }).map(({ x1, x2, y1, y2}) => ((x2 - x1) * (y1 + y2))).reduce((a, b) => a + b);
    
    if (sum > -1) return true;
    if (sum < 0) return false;
}

Take all the vertices as an array like this;

const vertices = [{x: 5, y: 0}, {x: 6, y: 4}, {x: 4, y: 5}, {x: 1, y: 5}, {x: 1, y: 0}];
isClockwise(vertices);

Solution for R to determine direction and reverse if clockwise (found it necessary for owin objects):

coords <- cbind(x = c(5,6,4,1,1),y = c(0,4,5,5,0))
a <- numeric()
for (i in 1:dim(coords)[1]){
  #print(i)
  q <- i + 1
  if (i == (dim(coords)[1])) q <- 1
  out <- ((coords[q,1]) - (coords[i,1])) * ((coords[q,2]) + (coords[i,2]))
  a[q] <- out
  rm(q,out)
} #end i loop

rm(i)
    
a <- sum(a) #-ve is anti-clockwise

b <- cbind(x = rev(coords[,1]), y = rev(coords[,2]))

if (a>0) coords <- b #reverses coords if polygon not traced in anti-clockwise direction

While these answers are correct, they are more mathematically intense than necessary. Assume map coordinates, where the most north point is the highest point on the map. Find the most north point, and if 2 points tie, it is the most north then the most east (this is the point that lhf uses in his answer). In your points,

point[0] = (5,0)

point[1] = (6,4)

point[2] = (4,5)

point[3] = (1,5)

point[4] = (1,0)

If we assume that P2 is the most north then east point either the previous or next point determine clockwise, CW, or CCW. Since the most north point is on the north face, if the P1 (previous) to P2 moves east the direction is CW. In this case, it moves west, so the direction is CCW as the accepted answer says. If the previous point has no horizontal movement, then the same system applies to the next point, P3. If P3 is west of P2, it is, then the movement is CCW. If the P2 to P3 movement is east, it's west in this case, the movement is CW. Assume that nte, P2 in your data, is the most north then east point and the prv is the previous point, P1 in your data, and nxt is the next point, P3 in your data, and [0] is horizontal or east/west where west is less than east, and [1] is vertical.

if (nte[0] >= prv[0] && nxt[0] >= nte[0]) return(CW);
if (nte[0] <= prv[0] && nxt[0] <= nte[0]) return(CCW);
// Okay, it's not easy-peasy, so now, do the math
if (nte[0] * nxt[1] - nte[1] * nxt[0] - prv[0] * (nxt[1] - crt[1]) + prv[1] * (nxt[0] - nte[0]) >= 0) return(CCW); // For quadrant 3 return(CW)
return(CW) // For quadrant 3 return (CCW)

Here's a simple Python 3 implementation based on this answer (which, in turn, is based on the solution proposed in the accepted answer)

def is_clockwise(points):
    # points is your list (or array) of 2d points.
    assert len(points) > 0
    s = 0.0
    for p1, p2 in zip(points, points[1:] + [points[0]]):
        s += (p2[0] - p1[0]) * (p2[1] + p1[1])
    return s > 0.0

find the center of mass of these points.

suppose there are lines from this point to your points.

find the angle between two lines for line0 line1

than do it for line1 and line2

...

...

if this angle is monotonically increasing than it is counterclockwise ,

else if monotonically decreasing it is clockwise

else (it is not monotonical)

you cant decide, so it is not wise