How to discover winding direction of Polygons?

[TJ1]

How do you programmatically detect the winding direction (order of coordinates - clockwise (CW) or counter-clockwise (CCW)) of a series of coordinates of a closed Polygon?

The winding direction of a Polygon determines which face is front and which back. Front face is indicated by a CCW winding direction.

This is absorbing and annoying me that something looking so easy to the human eye, is incredibly complicated to decide mathematically

I doubt I'll need it now I know the issue as it relates to GE, but its bugging me because I'm one of those kind of people that refuses to accept that these things have to be complicated - if I experiment enough and use my imagination I can usually come up with a sneaky way of getting a result without too much hassle.

If you know how its done 'officially' please tell me.

I seem to have come up with a way to determine it without getting into complicated maths with vectors and angles, but before I open my mouth and put my foot in it to prove how idiotic I am I want some ideas of the mathematical side of the problem

[Valery35]

On a sign on the area

A=1/2*sum((Xk+Xk+1)*(Yk-Yk+1)),
k=1...n,
Xk+1=X0, Yk+1=Y0

http://algolist.manual.ru/maths/geom/polygon/gif/polsquare.gif

[TJ1]

Yes, but you try writing fast code to calculate the area of a complex Polygon with, say, 2000 coordinates!!

Imagine you've got several hundred Polygons to determine the winding direction for... you're requiring some serious number crunching.

[TJ1]

Here's my current algorithm:

1. Enter the Points in an ordered array, with index 0 being the first in the list, and index n being the last.
2. Scan all the points to discover the four extents: most northly, easterly, southerly, westerly.
3. Record the indexes of the Points at the four extents (N, S, E, W).
4. Renumber the array so that the most northerly Point is at index 0 and the rest follow in original order, keeping track of the new indexes of N, S, E, W.
5. Compare the new index of the most Westerly Point with the new index of the most Easterly Point.
6. If the Westerly index is less than the Easterly index, the Polygon is counter-clockwise (CCW), otherwise if the Easterly index is less than the Westerly index, the Polygon is clockwise (CW).

I guess I need to generate some very complex models to see if this falls down, and if so, when.

[Valery35]

..4. Renumber the array so that the most northerly Point is at index 0 and the rest follow in original order, keeping track of the new indexes of N, S, E, W...

imho
If you use sorting it is algorithm of type (n^2)

Classical systems all the same count, as I have resulted.
They sometimes lead numbers to the integer (ESRI ArcGIS).

[TJ1]

The thing is you don't have to physically re-order the array, and you're not going to sort it.

You iterate it once, checking for a larger North South East and West than you've already found, and updating those if the current point exceeds them.
Code:

---

$maxNorth = $maxEast = $maxSouth = $maxWest = 0;
$maxNorthIdx = $maxEastIdx = $maxSouthIdx = $maxWestIdx = 0;
for($i=0; $i<count($points); $i++) {

 if ($points[$i]['longitude'] < $maxWest) {
   $maxWest = $points[$i]['longitude'];
   $maxWestIdx = $i;
 }
 else if ($points[$i]['longitude'] > $maxEast) {
   $maxEast = $points[$i]['longitude'];
   $maxEastIdx = $i;
 }

 if ($points[$i]['latitude'] < $maxSouth) {
   $maxSouth = $points[$i]['latitude'];
   $maxSouthIdx = $i;
 } 
 else if ($points[$i]['latitude'] > $maxNorth) {
   $maxNorth = $points[$i]['latitude'];
   $maxNorthIdx = $i;
 }
}
$offset = $maxNorthIdx; // alias all indexes by this amount

---

Then you alias the starting point from index '0' to Point North - so for example, if North is originally at index 4500, and the total number of points in the array is 5000, all point indexes simply have an integer offset subtracted from them, with simple logic to detect the wrap around using modulo array size:

0 Point A
1 Point B
2 Point C
...
500 Point West
...
2000 Point South
...
3000 Point East
...
4500 Point North
...
5000 Point Z

The offset is therefore 4500

So for any given Point, it's index is:

$index = ($index - $offset) < 0 ? ($index - $offset) + 5000 : $index - $offset;

So with no re-arranging of the array, we can calculate the West and East indexes and therefore the winding direction:

Code:

---

$West = ($maxWestIdx-offset) < 0 ? ($maxWestIdx-$offset)+5000 : $maxWestIdx-$offset;
$East = ($maxEastIdx-offset) < 0 ? ($maxEastIdx-$offset)+5000 : $maxEastIdx-$offset;

$CCW = $West < $East ? 1 : 0; // set boolean  flags to indicate winding direction
$CW = $West < $East ? 0 : 1;

---

This all supposes the algorithm is robust. For the cases I'm working with it works and is very fast, but I'm sure someone out there can show it falls down in some cases - the thing is, finding out when that happens, and whether those cases can be predicted and corrected.

As long as all the points are on or near the same plane it seems as it it will be reliable enough - certainly for working with GE Polygons. And it meets my criteria - easy to implement and fast.

As to the more general mathematical proofs, I'll leave that for the people that find that stuff entertaining

[Valery35]

a=0
do i=1 to count-1
a=a+(x(i)+x(i+1))*(y(i)-y(i+1))
enddo
if a>0 then else


??? =)

[TJ1]

Mine beats yours

Yours needs floating-point multiplication - mine only needs floating-point comparison and copy

Can you tell I'm a machine-code programmer at heart - counting the CPU clock cycles

[Valery35]

Aha. But do not ask then how to count up weights for a transparency

[TJ1]

Hey Valery... thanks for this algorithm.

I've used it in boolean KMLDocument->isWindingDirectionCCW(KMLElement coordinates) rather than trust my 'faster' method and it seems to perform well