331

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).

Machavity
  • 30,841
  • 27
  • 92
  • 100
Stécy
  • 11,951
  • 16
  • 64
  • 89
  • 19
    PLEASE NOTE: The accepted answer, and many answers after it, require a lot of additions and multiplications (they are based on area calculations that end negative or positive; e.g. "shoelace formula"). Before implementing one of those, consider [lhf's answer](https://stackoverflow.com/a/1180256/199364), which is simpler/quicker - based on [wiki - orientation of simple polygon](https://en.wikipedia.org/wiki/Curve_orientation#Orientation_of_a_simple_polygon). – ToolmakerSteve Feb 05 '19 at 20:51
  • 3
    I always think of it in terms of the cross product of two adjacent vectors. If I walk around the perimeter of the polygon my head points out of the plane. I cross the out of plane vector into my walking direction vector to get the third direction in my coordinate system. If that vector points so that the interior is on my left it's counterclockwise; if the interior is on my right it's clockwise. – duffymo Oct 04 '19 at 12:18
  • Without in any way diminishing the interest of the question; for those coming here and who don't want to reinvent the wheel, you can checkout the [GEOS/libgeos](https://libgeos.org/) source code as explained [hereunder](https://stackoverflow.com/a/75120030/6630397) along with a small Python snippet using [Shapely](https://shapely.readthedocs.io), where this orientation check is available. – swiss_knight Jan 14 '23 at 17:48

25 Answers25

503

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, (x2 − x1)(y2 + y1). 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
Roberto Bonvallet
  • 31,943
  • 5
  • 40
  • 57
Beta
  • 96,650
  • 16
  • 149
  • 150
  • 37
    It's calculus applied to a simple case. (I don't have the skill to post graphics.) The area under a line segment equals its average height (y2+y1)/2 times its horizontal length (x2-x1). Notice the sign convention in x. Try this with some triangles and you'll soon see how it works. – Beta Jul 27 '09 at 14:20
  • I'd suggest proving it first for a triangle, then showing that any polygon can be split into triangles in a way that doesn't change its "clockwiseness". – Beta Aug 31 '09 at 15:44
  • Is there a particular name to this method? – Mark Ingram Nov 09 '09 at 17:19
  • 1
    Now that I think of it, the method of calculating area using a counterclockwise curve can be derived from Green's Theorem. I saw the method (without derivation) years ago-- I may have come up with idea of using it as a clockwiseness test myself, but it's a pretty obvious jump. – Beta Nov 09 '09 at 19:17
  • 1
    Is there a similar elegant solution for polygons in 3D space? – the Ritz Mar 12 '13 at 13:47
  • 1
    Found it here: http://stackoverflow.com/questions/12642256/python-find-area-of-polygon-from-xyz-coordinates – the Ritz Mar 12 '13 at 13:54
  • 1
    It also generalizes nicely to 3-D polytopes: for each face compute the signed area in x-y plane as described in this answer, multiply by the mean value of z, and sum it up for all faces. The absolute value of the result will be the polytope's volume. – Michael Aug 27 '13 at 20:39
  • 93
    A minor caveat: this answer assumes a normal Cartesian coordinate system. The reason that's worth mentioning is that some common contexts, like HTML5 canvas, use an inverted Y-axis. Then the rule has to be flipped: if the area is *negative*, the curve is clockwise. – LarsH Oct 11 '13 at 20:49
  • 2
    This method is wrong! I implemented it ad is works in 98% cases. But not always. [What worked for me (in every case) is this link](http://en.wikipedia.org/wiki/Curve_orientation). – Mr.Qbs Dec 13 '14 at 16:49
  • 2
    @Mr.Qbs: Please give us the simplest counterexample you know. – Beta Dec 14 '14 at 00:21
  • I checked my implementation to show you that it's properly done and to give you counterexample. Then I saw your coment "I'd suggest proving it for a triangle" so I tried it and it works! When i have 3 points (2 edges of my big polygon) I need to calculate 3 edges to get properly result (ab, bc, ca - imaginary edge). But when I try it only for 2 "real" edges (that realy belongs to my polygon) result sometimes is bad (its hard to give example - i compute 3d objects and I can see the bad arc quickly on modified model - sometimes arc is made in bad direction). – Mr.Qbs Dec 14 '14 at 11:38
  • 16
    @Mr.Qbs: So my method works, but if you *skip a vital part*, then it doesn't work. This is not news. – Beta Dec 14 '14 at 14:03
  • 12
    @Mr.Qbs: You always have to link the last point to the first one. If you have N points numbered from 0 to N-1, then you must calculate: `Sum( (x[(i+1) mod N] - x[i]) * (y[i] + y[(i+1) mod N]) )` for i = 0 to N-1. I.e., must must take the index Modulo N (`N ≡ 0`) The formula works only for **closed** polygons. Polygons have no imaginary edges. – Olivier Jacot-Descombes Jul 20 '15 at 14:45
  • 1
    Can the result be zero? And if it can, then what orientation is it? (It looks like some special cease). – Dan M. Aug 02 '16 at 17:58
  • 2
    @DanM. If the result is zero, that means that the positive and negative areas cancel out, as in a figure-eight. – Beta Aug 03 '16 at 02:23
  • 1
    Btw, It seems you only need to calculate sum of `xi + 1*yi - xi*yi + 1` to determine the answer. – Dan M. Aug 13 '16 at 22:30
  • 1
    @DanM.: True, but notice that my formula has only one multiplication, while yours has two. – Beta Aug 14 '16 at 19:54
  • I'm trying this out with a clockwise diamond, and getting a negative (counter-clockwise) sum result, and a positive sum for counter-clockwise. CLOCKWISE DIAMOND: `point[0] = 0,5 (5-0)(0+5) = 25, point[1] = 5,0 (10-5)(5+0) = 25, point[2] = 10,5 (5-10)(5+10) = -75, point[3] = 5,10 (0-5)(10+5) = -75` == -100. COUNTER DIAMOND: `point[0] = 0,5 (5-0)(5+10) = 75, point[1] = 5,10 (10-5)(10+5) = 75, point[2] = 10,5 (5-10)(5+0) = -25, point[3] = 5,0 (0-5)(5+0) = -25` == +100 Is this a special case? – Matt Fiocca Jan 12 '17 at 02:58
  • @MattFiocca: You are using the words "clockwise" and "counter-clockwise" in the opposite sense of everyone else. – Beta Jan 12 '17 at 03:09
  • @Beta: just so i'm not going crazy, is everyone else using clockwise as counter-clockface? Here's a link to my diagram with calculations. https://s3-us-west-2.amazonaws.com/mattfiocca/polydirection.png. – Matt Fiocca Jan 12 '17 at 05:28
  • @MattFiocca: The convention in mathematics is that the Y coordinate increases upward; your Y coordinates increase downward. – Beta Jan 12 '17 at 23:41
  • @Beta: ah thanks! I knew I was missing something vital, having come from the backward world of web and apps. – Matt Fiocca Jan 13 '17 at 03:33
  • This works only with all positive co-ordinates. I tried with Polygons in X,-Z Plane (all Z values are negative in this case.). Only way I could get this logic to work in this case is by translating the polygons to X,+Z plane. i.e, by offsetting the Z values of all co-ordinates with a constant such that all values are positive. Please correct me if am wrong or if there is any simpler alternative. – Sreeraj Mar 17 '17 at 07:03
  • @Sreeraj: You are mistaken, but a comment thread is not the place to discuss it. If I have time I'll make a chat room. – Beta Mar 18 '17 at 22:08
  • 12
    This http://blog.element84.com/polygon-winding.html explains in simple english why this solution works. – David Zorychta May 17 '17 at 01:16
  • 1
    Alternatively it can be calculated as sum of `(point.x + nextPoint.x) * (point.y - nextPoint.y)`. – RunninglVlan Nov 26 '17 at 16:14
  • 1
    Would this work for lat/lon coordinates assuming you're not near/crossing a pole, and all polygons that cross 180 are defined with numbers greater than 180? – Sandy Gifford Apr 16 '18 at 11:36
  • 1
    @Macmee: the link given by you is broken, archive: https://web.archive.org/web/20171212175910/blog.element84.com/polygon-winding.html – Ans Nov 27 '18 at 14:51
  • 1
    Ah, with a clever approach that cancels the extra XnYn terms over the sum, this is a cross product. Sweet! – Brian Carcich Dec 04 '18 at 16:57
  • @LarsH unless it's Y coordinate is inverted resulting in a negative Y value, wouldn't positive sum of edges still means the polygon is clockwise? – Chanwoo Park Apr 23 '20 at 05:50
  • 1
    @ChanwooPark: Good question. I can see why that might appear to be the case, since `y1 + y2` wouldn't change sign just from inverting the direction of the y axis. But I don't think so. Try it for a simple triangle. (I don't have time to try it now.) – LarsH Apr 23 '20 at 14:21
  • 1
    Ok I tried it with the example given in the answer. If you flip the y-axis to go from 5 to 0 instead of 0 to 5, you end up with 45 instead of -44. So no, a positive sum doesn't mean the polygon is clockwise, when the y-axis is flipped. – LarsH Apr 23 '20 at 14:28
  • 1
    @ChanwooPark Think of how the algorithm works: if edges that go leftward (x2 < x1) have higher y-values overall than edges that go rightward (x2 > x1), the negative terms will outweigh the positive, and in a Cartesian coordinate system this means the polygon goes counterclockwise. – LarsH Apr 23 '20 at 15:28
  • 1
    Whereas in an inverted y-coordinate system, edges that would have had higher y-values than the others now have lower y-values. So for a counterclockwise polygon, the leftward edges now have lower y-values, and the positive terms of the rightward edges will outweigh them, producing a positive result. – LarsH Apr 23 '20 at 15:31
  • @ChanwooPark Thx for the question. I hadn't really understood the algorithm until you provoked me to think it thru. – LarsH Apr 24 '20 at 15:24
  • Can this be generalized to 3 dimensions, x, y, z? – stenfeio Feb 15 '21 at 22:59
  • 1
    @stenfeio: Yes, but you must think carefully about what "clockwise" means in 3D. If you don't see how, you can post ta Question; it's too much to add to this one. – Beta Feb 16 '21 at 01:09
  • @Beta, posted here: https://stackoverflow.com/questions/66226478/how-do-you-determine-if-a-list-of-of-points-in-3d-space-are-in-clock-wise-order – stenfeio Feb 16 '21 at 14:41
88

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:

Michael
  • 780
  • 6
  • 14
lhf
  • 70,581
  • 9
  • 108
  • 149
  • 11
    This is also explained in http://en.wikipedia.org/wiki/Curve_orientation. The point is that the the found point must be on the convex hull, and it's only necessary to look locally at a single point on the convex hull (and its immediate neighbors) to determine the orientation of the whole polygon. – M Katz Jun 02 '15 at 06:30
  • 2
    Shocked and awed this hasn't received more upvotes. For simple polygons (*which is most polygons in some fields*), this answer yields an `O(1)` solution. All other answers yield `O(n)` solutions for `n` the number of polygon points. For even deeper optimizations, see the [*Practical Considerations*](https://en.wikipedia.org/wiki/Curve_orientation#Practical_considerations) subsection of Wikipedia's fantastic [Curve orientation](http://en.wikipedia.org/wiki/Curve_orientation) article. – Cecil Curry Jun 03 '17 at 07:39
  • 16
    **_Clarification:_** this solution is `O(1)` only if either **(A)** this polygon is convex (in which case any arbitrary vertex resides on the convex hull and hence suffices) *or* **(B)** you already know the vertex with the smallest Y coordinate. If this is *not* the case (i.e., this polygon is non-convex and you don't know anything about it), an `O(n)` search is required. Since no summation is required, however, this is still dramatically faster than any other solution for simple polygons. – Cecil Curry Jun 03 '17 at 07:49
  • An implementation of this answer: [c# code to find corner vertex and calculate determinant of angle at that vertex.](https://stackoverflow.com/a/54546353/199364) – ToolmakerSteve Feb 06 '19 at 03:51
  • 3
    @CecilCurry I think your 2nd comment explains why this hasn't received more upvotes. It yields wrong answers in certain scenarios, without any mention of those limitations. – LarsH May 02 '20 at 15:16
61

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 * (x1*y2 - x2*y1 + x2*y3 - x3*y2 + ... + xn*y1 - x1*yn)

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

StayOnTarget
  • 11,743
  • 10
  • 52
  • 81
Sean the Bean
  • 5,222
  • 5
  • 38
  • 40
  • Was that a typo in your signed area formula above? It ends with "xn*y1 - x1*yn"; when I believe it should be "x_n y_{n+1} - y_n x_{n-1}" (in LaTeX, at least). On the other hand, it's been ten years since I took any linear algebra classes. – Michael Macha May 25 '15 at 20:25
  • Nope. If you check the [source](http://mathworld.wolfram.com/PolygonArea.html), you'll see that the formula does in fact reference the first point again in the last term (y1 and x1). (Sorry, I'm not very familiar with LaTeX, but I formatted the subscripts to make them more readable.) – Sean the Bean May 26 '15 at 17:38
  • I used this solution and it worked perfectly for my use. Note that if you can plan ahead and spare and extra two vectors in your array, you can get rid of the comparison (or %) by adding the first vector at the tail of the array. That way you simply loop over all the elements, except the last one (length-2 instead of length-1). – Eric Fortier Jan 16 '19 at 02:03
  • 3
    @EricFortier - FWIW, rather than resize a possibly large array, an efficient alternative is for each iteration to save its point as `previousPoint` for next iteration. Before starting loop, set `previousPoint` to array's last point. Trade off is extra local variable copy but fewer array accesses. And most importantly, don't have to touch the input array. – ToolmakerSteve Feb 05 '19 at 20:16
  • 3
    @MichaelEricOberlin - its necessary to *close* the polygon, by including the line segment from last point to first point. (A correct calculation will be the same, no matter which point starts the closed polygon.) – ToolmakerSteve Feb 05 '19 at 20:20
  • @ToolmakerSteve - Good point. Unless you know the size beforehand and can create the array with the extra element, your solution is fast and prevents resizing. – Eric Fortier Feb 10 '19 at 16:45
51

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.

Charles Bretana
  • 143,358
  • 22
  • 150
  • 216
  • 3
    Actually, this solution is a different solution than the accepted solution. Whether they are equivalent or not is a question I am investigating, but I suspect they are not... The accepted answer calculates the area of the polygon, by taking the difference between the area under the top edge of the polygon and the area under the bottom edge of the polygon. One will be negative (the one where you are traversing from left to right), and the other will be negative. When traversing clockwise, The upper edge is traversed left to right and is larger, so the total is positive. – Charles Bretana Apr 06 '13 at 14:29
  • 1
    My solution measures the sum of sines of the changes in edge angles at each vertex. This will be positive when traversing clockwise and negative when traversing counter clockwise. – Charles Bretana Apr 06 '13 at 14:30
  • 2
    It seems with this approach you DO need to take the arcsin, unless you assume convexity (in which case you need only check one vertex) – agentp May 30 '13 at 16:54
  • No, you do not need to take the arcsine. You do need to divide by the magnitudes of the two vector, to eliminate scaling, (where one very large pair od segments that bends one way a small angle contributes more tot the sum than many smaller segments that bend the other way by a larger angle. But the arcsine just converts one measure of an angle to another. – Charles Bretana May 05 '15 at 23:53
  • Why all the expensive math functions? Wouldn't it be faster to only calculate the z-part of a 3D cross product by doing a1 * b2 - a2 * b1 (as the x- and y-parts will be zero anyway) and checking the sign of the result? No trigonometric functions are needed and no sqrt. – marsbear Jan 19 '18 at 20:24
  • @marsbear, That will only work if all vectors are the same magnitude (or, equivalently, have first been normalized). If they are of different lengths, then you will get the wrong answer., because a clockwise angle between two long segments of the polygon will contribute much more to the overall sum than a counterclockwise angle between two very small segments. – Charles Bretana Jan 19 '18 at 21:51
  • But shouldn't it be enough to "add up" the signs of the cross products? – marsbear Jan 24 '18 at 10:25
  • No, all that would do is "count" the number of turns to the right minus the number of turns to the left. If you start pointing East, and take ninety one degree turns to the right, (totaling 90 degrees of right), each an inch long, you will be pointing south, and then take five 90 degree turns to the left, you will be pointing east again, and with appropriate adjustment of the length of each segment, (make the lefts longer) reconnect to the start point. The overall closed curve is clearly counterclockwise, but just counting the signs would count 90 to the right and only 5 to the left. – Charles Bretana Jan 24 '18 at 15:52
  • I have reformulated your answer (to make it more precise mathematically) until your sentence "The magnitude of this value is a measure of the sine of the angle between the 2 original vectors". You should edit that sentence and the remaining ones to take into account my changes. The reason I didn't edit that sentence and the remaining ones is that I am not fully understanding what you mean those. – nbro Feb 24 '18 at 23:46
  • 3
    You DO need to take the arcsin. Try it on a bunch of random non-convex polygons, and you will find the test will fail for some polygons if you don't take the arcsin. – Luke Hutchison Jun 08 '18 at 09:07
  • No, You don't. Whatever test you are running, you are running it wrong. The value before taking the arcsine IS the sin of the original angle of rotation, converting it from a value between -1 and +1 to an actual angle between -Pi and +Pi does nothing to change the answer. And as long as you divide by the product of the magnitudes of the two vectors, the magnitudes will be normalized, and their sum will indicate the overall direction of curvature. – Charles Bretana Jun 12 '18 at 14:48
  • 1
    @CharlesBretana - while I have not run Luke's test, I believe he is correct. That is the nature of *summing* combined with a *nonlinear* scale [without arcsin vs. with arcsin]. Consider what marsbear suggested, that you correctly rejected. He suggested that you "just count", and you pointed out that a handful of large values could outweigh a large number of small values. Now consider arcsin of each value vs not. Isn't it still the case that failing to take arcsin gives incorrect weight to each value, therefore has same flaw (though much less so)? – ToolmakerSteve Feb 05 '19 at 20:03
36

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.


Test (See discussion with @WDUK)

public static bool IsClockwise(IList<(double X, double Y)> vertices)
{
    double sum = 0.0;
    var v1 = vertices[^1];
    for (int i = 0; i < vertices.Count; i++) {
        var v2 = vertices[i];
        sum += (v2.X - v1.X) * (v2.Y + v1.Y);
        Console.WriteLine($"(({v2.X,2}) - ({v1.X,2})) * (({v2.Y,2}) + ({v1.Y,2})) = {(v2.X - v1.X) * (v2.Y + v1.Y)}");
        v1 = v2;
    }
    Console.WriteLine(sum);
    return sum > 0.0;
}

public static void Test()
{
    Console.WriteLine(IsClockwise(new[] { (-5.0, -5.0), (-5.0, 5.0), (5.0, 5.0), (5.0, -5.0) }));

    // infinity Symbol
    //Console.WriteLine(IsClockwise(new[] { (-5.0, -5.0), (-5.0, 5.0), (5.0, -5.0), (5.0, 5.0) }));
}
Olivier Jacot-Descombes
  • 104,806
  • 13
  • 138
  • 188
  • 2
    You can avoid the costly `%` and avoid branching too by setting `v1 = vertices[vertices.Count-1]` before the loop starts, use `v2 = vertices[i];` then after the addition to the `sum` do `v1 = v2`. – Michel Rouzic Feb 18 '22 at 15:09
  • @MichelRouzic how is `(i+1) % (count-1)` a branch isn't that just a math operation for the next index? Why would this involve branching ? – WDUK Apr 28 '22 at 08:24
  • It took me some time to remember what I meant by this, it's not branching, I meant that if anyone tried to do something like `(i+1) % (count-1)` without `%` they'd go for something like `i+1 < count ? i+1 : 0` like [this answer](https://stackoverflow.com/a/49035316/1675589) for instance, when there's a way that avoids both modulo and branching. – Michel Rouzic Apr 29 '22 at 12:37
  • This doesn't work if the points are in symmetrical positions Example of 4 points for a square: (-5,-5)(-5,5)(5,5)(5,-5) you will get zero when they are clearly clockwise ordered... – WDUK Dec 21 '22 at 09:15
  • @WDUK, I get a sum of 200: (-5 - 5) * (-5 + -5) = 100, (-5 - -5) * (5 + -5) = 0, (5 - -5) * (5 + 5) = 100, (5 - 5) * (-5 + 5) = 0. The sum is 0 if you invert the two last points, in which case get a shape looking like an infinity symbol. – Olivier Jacot-Descombes Dec 21 '22 at 11:10
  • `(-5 - 5) * (-5 + -5) = 100` ? The first two points are (-5,-5) and (-5,5) So v2.x - v1.x will always be 0 so how did you get 100 for the first two points? `((-5) - (-5)) * ((-5) - 5) = 0` – WDUK Dec 21 '22 at 22:04
  • @WDUK, I added the test to my answer. Note that the first sum is calculated from the last (`v1 = vertices[^1]`) and first (`v2`) point. So the zero you mention is printed as the second line in the test. – Olivier Jacot-Descombes Dec 22 '22 at 13:05
7

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:

mpen
  • 272,448
  • 266
  • 850
  • 1,236
6

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.

Steve Gilham
  • 11,237
  • 3
  • 31
  • 37
  • Seeing as how the cross product basically boils down to a positive scaling factor times the sine of the angle, it's probably better to just do a cross product. It'll be faster and less complicated. – ReaperUnreal Jul 22 '09 at 14:44
4

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;
  }
})();
Steve Jansen
  • 9,398
  • 2
  • 29
  • 34
  • Trying this I get exactly the opposite result, a polygon drawn in clockwise order yields a negative area, while one drawn counter clockwise yields positive. In either case, this snippet is still super useful 5yrs on, thank you. – Cameron Roberts Sep 27 '18 at 20:16
  • @CameronRoberts The norm (see IETF in particular for geoJson) is to follow the 'right-hand rule'. I guess that Google is complaining with. In that case outer ring must be counterclockwise (yielding positive area), and inner rings (holes) are winding clockwise (negative area to be removed from main area). – allez l'OM Dec 03 '19 at 13:46
4

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);
}
nbro
  • 15,395
  • 32
  • 113
  • 196
MSS
  • 3,520
  • 24
  • 29
  • Openlayers is javascript based map management library like googlemaps and it's wrote and used in openlayers 2. – MSS Feb 26 '18 at 07:29
  • Can you explain a little bit what your code does, and why you're doing it? – nbro Feb 26 '18 at 10:25
  • @nbro this code implements the [lhf answer](https://stackoverflow.com/questions/1165647/how-to-determine-if-a-list-of-polygon-points-are-in-clockwise-order/1180256#1180256). It is easy to keep the non OpenLayer part in a pure javascript function by having _vertices_ directly as parameter. It works well, and could be adapted to the case of _multiPolygon_. – allez l'OM Dec 03 '19 at 14:31
4

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;
}
ToolmakerSteve
  • 18,547
  • 14
  • 94
  • 196
  • 2
    This appears to be for down-is-positive Y coordinates. Flip CW/CCW for standard coordinates. – Warwick Allison Feb 09 '19 at 00:37
  • `down-is-positive Y` *is* standard coordinates. – Tatarize May 08 '23 at 22:51
  • @Tatarize - that depends on the context. If no context is stated, it is generally assumed to be up. OTOH, since this is a computer graphics discussion, it could be argued either way. [Cartesian coordinate system](https://en.wikipedia.org/wiki/Cartesian_coordinate_system): *"In mathematics, physics, and engineering, the first axis is usually defined or depicted as horizontal and oriented to the right, and the second axis is vertical and oriented **upwards** [emphasis added]. (However, in some computer graphics contexts, the ordinate axis may be oriented downwards.)"* – ToolmakerSteve May 09 '23 at 00:43
  • It does depend a bit on context. Most graphing stuff including math graphs are `up-is-positive-y` but graphic systems which is a lot of computer science is are `down-is-positive-y` systems. There really isn't a "standard" there. And, in fact, in a lot of projects I make a point to write a coordinate system section within the documentation. – Tatarize May 09 '23 at 03:45
  • Exactly. Now we are in agreement. – ToolmakerSteve May 09 '23 at 17:46
2

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

nbro
  • 15,395
  • 32
  • 113
  • 196
Frederick
  • 1,271
  • 1
  • 10
  • 29
1

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

enter image description here

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.

nbro
  • 15,395
  • 32
  • 113
  • 196
Ian
  • 3,500
  • 1
  • 24
  • 25
  • @tibetty Can you explain why this method wouldn't work in many situations if the polygon is concave? – nbro Feb 25 '18 at 11:34
  • 1
    Please look at the last table in the wiki item reference in your post. It's easy for me to give a false example but hard to prove it. – tibetty Feb 28 '18 at 01:20
  • 1
    Please look at the last table in the wiki item reference in your post. It's easy for me to give a false example but hard to prove it. – tibetty Feb 28 '18 at 01:22
  • 1
    @tibetty is correct. You can't simply take any three points along the polygon; you might be in either a convex or concave region of that polygon. Reading wiki carefully, one must take three points *along the convex hull that encloses the polygon*. From "practical considerations": *"One does not need to construct the convex hull of a polygon to find a suitable vertex. A common choice is the vertex of the polygon with the smallest X-coordinate. If there are several of them, the one with the smallest Y-coordinate is picked. It is guaranteed to be [a] vertex of the convex hull of the polygon."* – ToolmakerSteve Feb 05 '19 at 21:07
  • 2
    Hence [lhf's earlier answer](https://stackoverflow.com/a/1180256/199364), which is similar, and references the same wiki article, but specifies such a point. [Apparently it doesn't matter whether one takes smallest or largest, x or y, as long as one avoids being in the middle; effectively one is working from one edge of the bounding box around the polygon, to guarantee in a concave region.] – ToolmakerSteve Feb 05 '19 at 21:11
1

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
nbro
  • 15,395
  • 32
  • 113
  • 196
0

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
nbro
  • 15,395
  • 32
  • 113
  • 196
Gianni Spear
  • 7,033
  • 22
  • 82
  • 131
0

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.

daniel
  • 157
  • 1
  • 1
  • 1
    True, but you misunderstand the concept of a polygon's winding order (clockwise or counter-clockwise). In an entirely convex polygon, the angle at all points will be clockwise or all will be counter-clockwise [as in your first sentence]. In a polygon with concave region(s), the "caves" will be in the opposite direction, but the polygon as a whole still has a well-defined interior, and is considered clockwise or counter-clockwise accordingly. See https://en.wikipedia.org/wiki/Curve_orientation#Orientation_of_a_simple_polygon – ToolmakerSteve Feb 05 '19 at 21:40
0

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();



        } 
    }
}
bradgonesurfing
  • 30,949
  • 17
  • 114
  • 217
0

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.

nbro
  • 15,395
  • 32
  • 113
  • 196
Marjan Moderc
  • 2,747
  • 23
  • 44
0

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.

nbro
  • 15,395
  • 32
  • 113
  • 196
  • *Maybe* this works if the polygon is entirely convex. It definitely is not reliable if there are any concave regions - its easy to pick a point that is on the "wrong" side of one of the edges of the cave, then connect it to that edge. Will get wrong answer. – ToolmakerSteve Feb 05 '19 at 21:22
  • It works even if the polygon is concave. The point needs to be inside that concave polygon. However I am not sure about complex polygon (did not test.) – Venkata Goli Feb 06 '19 at 22:28
  • "It works even if the polygon is concave." - Counterexample: poly (0,0), (1,1), (0,2), (2,2), (2,0). Line segment (1,1), (0, 2). If you pick an interior point within (1,1), (0,2), (1,2) to form triangle -> (1,1), (0,2), (0.5,1.5)), you get *opposite* winding than if you pick an interior point within (0,0), (1,1), (1,0) > (1,1), (0,2),(0.5,0.5). Those are both interior to the original polygon, yet have opposite windings. Therefore, one of them gives the wrong answer. – ToolmakerSteve Feb 07 '19 at 18:29
  • In general, if a polygon has any concave region, pick a segment in the concave region. Because it is concave, you can find two "interior" points that are on opposite sides of that line. Because they are on opposite sides of that line, the triangles formed have opposite windings. End of proof. – ToolmakerSteve Feb 07 '19 at 18:42
0

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
Toby Evetts
  • 123
  • 1
  • 9
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);
Ind
  • 41
  • 5
0

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
dez
  • 1
  • 1
0

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)
VectorVortec
  • 667
  • 7
  • 10
  • IMHO, it would be safer to stick to the fundamental math shown in [lhf's answer](https://stackoverflow.com/a/1180256/199364) - thank you for mentioning him. The challenge in reducing it to quadrants is that its a fair amount of work to *prove* that your formula is correct in all cases. Did you correctly calculate "more west"? In a concave polygon where *both* [1] and [3] are "west and south" of [2]? Did you correctly handle different lengths of [1] and [3] in that situation? I have no idea, whereas if I directly compute that angle (or its determinant), I'm using well-known formulas. – ToolmakerSteve Feb 05 '19 at 21:55
  • @ToolmakerSteve the if statements always work if the 3 points are convex. The if statements will return, then you'll get the right answer. The if statements will not return, if the shape is concave and extreme. That's when you have to do the math. Most images have one quadrant, so that part is easy. More than 99% of my subroutine calls are handled by the if statements. – VectorVortec Feb 06 '19 at 01:42
  • That doesn't address my concern. What is that formula? Is it the orientation determinant as given in the wiki link from lhf's answer? If so, then say so. Explain that what you are doing is doing quick checks that handle most cases, to avoid the standard math. If that is so, then your answer now makes sense to me. (Minor nit: would be easier to read if you used `.x` and `.y` of a struct, instead of `[0]` and `[1]`. I didn't know what your code was saying, first time I glanced at it.) – ToolmakerSteve Feb 06 '19 at 03:37
  • Since I did not have confidence in your approach, I [implemented lhf's approach](https://stackoverflow.com/a/54546353/199364); formula from his link. Slow part is *finding* appropriate vertex - O(N) search. Once found, the determinant is an O(1) operation, using 6 multiplies with 5 adds. That last part is what you've optimized; but you've done so by adding additional if-tests. I can't personally justify taking a non-standard approach - would need to verify each step is correct - But thank you for an interesting analysis of quadrants! – ToolmakerSteve Feb 06 '19 at 04:07
0

Javascript implementation of lhf's answer
(Again, this only works for simple polygons, e.i. not for figure 8's)

let polygon = [
  {x:5,y:0},
  {x:6,y:4},
  {x:4,y:5},
  {x:1,y:5},
  {x:1,y:0}
]
document.body.innerHTML += `Polygon ${polygon.map(p=>`(${p.x}, ${p.y})`).join(", ")} is clockwise? ${isPolygonClockwise(polygon)}`

let reversePolygon = []
polygon.forEach(point=>reversePolygon.unshift(point))
document.body.innerHTML += `<br/>Polygon ${reversePolygon.map(p=>`(${p.x}, ${p.y})`).join(", ")} is clockwise? ${isPolygonClockwise(reversePolygon)}`

function isPolygonClockwise (polygon) {
  // From http://www.faqs.org/faqs/graphics/algorithms-faq/ "How do I find the orientation of a simple polygon?"
  // THIS SOMETIMES FAILS if the polygon is a figure 8, or similar shape where it crosses over itself
  
  // Take the lowest point (break ties with the right-most). 
  if (polygon.length < 3) {
    return true // A single point or two points can't be clockwise/counterclockwise
  }
  let previousPoint = polygon[0]
  let lowestPoint = polygon[1]
  let nextPoint = polygon[2]
  polygon.forEach((point, index)=>{
    if (point.y > lowestPoint.y || (point.y === lowestPoint.y && point.x > lowestPoint.x)) { // larger y values are lower, in svgs
      // Break ties with furthest right
      previousPoint = polygon[(index-1) >= (0)                ? (index-1) : (polygon.length-1)]
      lowestPoint = polygon[index]
      nextPoint = polygon[(index+1)     <= (polygon.length-1) ? (index+1) : (0)]
    }
  })
  // Check the angle between the previous point, that point, and the next point.
  // If the angle is less than PI radians, the polygon is clockwise
  let angle = findAngle(previousPoint, lowestPoint, nextPoint)
  return angle < Math.PI
}

function findAngle(A,B,C) {
  var AB = Math.atan2(B.y-A.y, B.x-A.x);
  var BC = Math.atan2(C.y-B.y, C.x-B.x);
  if (AB < 0) AB += Math.PI*2
  if (BC < 0) BC += Math.PI*2
  return BC-AB;
}
Richard Henage
  • 1,758
  • 1
  • 3
  • 10
0

For those who don't want to "reinvent the wheel", I think it's worth mentioning that this check is implemented in a nice Python package called Shapely (github) (which is based on the GEOS C/C++ library):

Shapely is a BSD-licensed Python package for manipulation and analysis of planar geometric objects. It is using the widely deployed open-source geometry library GEOS (the engine of PostGIS, and a port of JTS). Shapely wraps GEOS geometries and operations to provide both a feature rich Geometry interface for singular (scalar) geometries and higher-performance NumPy ufuncs for operations using arrays of geometries. Shapely is not primarily focused on data serialization formats or coordinate systems, but can be readily integrated with packages that are.

Source: https://shapely.readthedocs.io/en/stable/

A small example given the OP coordinates:

import numpy as np
from shapely.geometry import Polygon

points = np.array([
    (5,0),
    (6,4),
    (4,5),
    (1,5),
    (1,0)
])

P = Polygon(points)

This is the newly constructed polygon:

import matplotlib.pyplot as plt

x,y = P.exterior.coords.xy
plt.plot(x,y)
plt.axis('equal')
plt.grid()
plt.show()

enter image description here

And you can directly use the is_ccw property of a LinearRing to check if the polygon is CW or CCW:

type(P.exterior)
>: shapely.geometry.polygon.LinearRing

P.exterior.is_ccw
>: True

If reversed:

points = np.flipud(points)
points
>: 
array([[1, 0],
       [1, 5],
       [4, 5],
       [6, 4],
       [5, 0]])


P1 = Polygon(points)

P1.exterior.is_ccw
>: True

Documentation and references for further reading:

swiss_knight
  • 5,787
  • 8
  • 50
  • 92
-4

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

ufukgun
  • 6,889
  • 8
  • 33
  • 55