Jak ustalić, czy lista punktów wielokąta jest zgodna z ruchem wskazówek zegara?

Mając listę punktów, jak znaleźć, czy są one w kolejności zgodnej z ruchem wskazówek zegara?

Na przykład:

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

powiedziałby, że jest przeciwny do ruchu wskazówek zegara (lub w przypadku niektórych osób przeciwnie do ruchu wskazówek zegara).

Niektóre z sugerowanych metod zawiodą w przypadku wielokąta niewypukłego, takiego jak półksiężyc. Oto prosty, który będzie działał z wielokątami niewypukłymi (zadziała nawet z wielokątem przecinającym się jak ósemka, informując, czy jest to w większości zgodnie z ruchem wskazówek zegara).

Suma ponad krawędziami, (x ₂ - x ₁ ) (y ₂ + y ₁ ). Jeśli wynik jest dodatni, krzywa jest zgodna z ruchem wskazówek zegara, a jeśli jest ujemna, krzywa jest przeciwna do ruchu wskazówek zegara. (Wynikiem jest dwukrotność zamkniętego obszaru z konwencją +/-).

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

Produkt przekroju mierzy stopień prostopadłej-ności dwóch wektorów. Wyobraź sobie, że każda krawędź twojego wielokąta jest wektorem w płaszczyźnie xy trójwymiarowej (3-D) przestrzeni xyz. Zatem iloczynem krzyżowym dwóch kolejnych krawędzi jest wektor w kierunku Z (dodatni kierunek Z, jeśli drugi segment jest zgodny z ruchem wskazówek zegara, minus kierunek Z, jeśli jest przeciwny do ruchu wskazówek zegara). Wielkość tego wektora jest proporcjonalna do sinusa kąta między dwiema oryginalnymi krawędziami, więc osiąga maksimum, gdy są one prostopadłe, i zwęża się, aby znikać, gdy krawędzie są współliniowe (równoległe).

Tak więc dla każdego wierzchołka (punktu) wielokąta oblicz wielkość iloczynu krzyżowego dwóch sąsiednich krawędzi:

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

Więc Wytwórnia krawędzie kolejno jak
edgeAto segment od point0do point1i
edgeBpomiędzy point1do point2
...
edgeEznajduje się pomiędzy point4i point0.

Zatem wierzchołek A ( point0) znajduje się między
edgeE[Od point4do point0]
edgeA[Od point0do `point1 '

Te dwie krawędzie same są wektorami, których współrzędne xiy można określić, odejmując współrzędne ich punktów początkowego i końcowego:

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

A iloczyn z tych dwóch przylegających krawędzi jest obliczana przy użyciu determinantę poniższej macierzy, która jest wykonana przez umieszczenie współrzędne dwóch wektorów poniżej symboli reprezentujących trzy oś współrzędnych ( i, j, i k). Trzecia (zerowa) współrzędna wartościowana istnieje, ponieważ koncepcja produktu krzyżowego jest konstrukcją 3-D, dlatego rozszerzamy te wektory 2-D na 3-D, aby zastosować produkt krzyżowy:

 i    j    k 
-4    0    0
 1    4    0    

Biorąc pod uwagę, że wszystkie produkty krzyżowe wytwarzają wektor prostopadły do ​​płaszczyzny dwóch wektorów, które są mnożone, wyznacznik powyższej macierzy ma tylko kskładnik (lub oś Z).
Wzór na obliczenie wielkości kskładnika w osi Z wynosi
a1*b2 - a2*b1 = -4* 4 - 0* 1 = -16

Wielkość tej wartości ( -16) jest miarą sinusoidy kąta między 2 oryginalnymi wektorami, pomnożonej przez iloczyn wielkości 2 wektorów.
W rzeczywistości inna formuła dla jego wartości to
A X B (Cross Product) = |A| * |B| * sin(AB).

Aby więc wrócić do miary kąta, należy podzielić tę wartość ( -16) przez iloczyn wielkości dwóch wektorów.

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

Więc miara grzechu (AB) = -16 / 16.4924=-.97014...

Jest to miara tego, czy następny segment po wierzchołku wygiął się w lewo lub w prawo i o ile. Nie ma potrzeby stosowania sinusoidy. Wszystko, na czym nam zależy, to jego wielkość i oczywiście jej znak (pozytywny lub negatywny)!

Zrób to dla każdego z pozostałych 4 punktów wokół zamkniętej ścieżki i dodaj wartości z tego obliczenia dla każdego wierzchołka.

Jeśli końcowa suma jest dodatnia, poszedłeś zgodnie z ruchem wskazówek zegara, ujemny, przeciwnie do ruchu wskazówek zegara.

Sądzę, że to dość stare pytanie, ale i tak zamierzam wyrzucić inne rozwiązanie, ponieważ jest proste i nie wymaga matematyki - używa tylko podstawowej algebry. Oblicz podpisany obszar wielokąta. Jeśli jest ujemny, punkty są w kolejności zgodnej z ruchem wskazówek zegara, a jeśli są dodatnie, są przeciwne do ruchu wskazówek zegara. (Jest to bardzo podobne do rozwiązania Beta.)

Oblicz podpisany obszar: A = 1/2 * (x ₁ * y ₂ - x ₂ * y ₁ + x ₂ * y ₃ - x ₃ * y ₂ + ... + x _n * y ₁ - x ₁ * y _n )

Lub w pseudokodzie:

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

Pamiętaj, że jeśli sprawdzasz tylko kolejność, nie musisz zawracać sobie głowy dzieleniem przez 2.

Źródła: http://mathworld.wolfram.com/PolygonArea.html

Znajdź wierzchołek o najmniejszym y (i największym x, jeśli istnieją powiązania). Niech wierzchołek będzie, Aa poprzedni wierzchołek na liście będzie, Ba następny wierzchołek na liście będzie C. Teraz obliczyć znak iloczynu krzyżowego ABi AC.

Bibliografia:

• Jak znaleźć orientację prostego wielokąta? w często zadawanych pytaniach: comp.graphics.algorytmy .

• Orientacja krzywej na Wikipedii.

Oto prosta implementacja algorytmu C # na podstawie tej odpowiedzi .

Załóżmy, że mamy Vectortyp mający Xi Ywłaściwości typu 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;
}

%jest operatorem modulo lub reszty wykonującym operację modulo, która ( według Wikipedii ) znajduje resztę po podzieleniu jednej liczby przez drugą.

Zacznij od jednego z wierzchołków i obliczyć kąt z każdej strony.

Pierwszy i ostatni będzie wynosił zero (więc pomiń te); dla reszty, sinus kąta będzie podany przez iloczyn krzyżowy normalizacji do długości jednostkowej (punkt [n] -punkt [0]) i (punkt [n-1] -punkt [0]).

Jeśli suma wartości jest dodatnia, wówczas wielokąt jest rysowany w kierunku przeciwnym do ruchu wskazówek zegara.

Dla tego, co jest warte, użyłem tego miksu do obliczenia kolejności nawijania dla aplikacji Google Maps API v3.

Kod wykorzystuje efekt uboczny obszarów wielokątów: kolejność nawijania wierzchołków w kierunku zgodnym z ruchem wskazówek zegara daje obszar dodatni, natomiast kolejność nawijania w kierunku przeciwnym do ruchu wskazówek zegara tych samych wierzchołków daje ten sam obszar jako wartość ujemną. Kod korzysta również z pewnego rodzaju prywatnego interfejsu API w bibliotece geometrii Map Google. Czułem się komfortowo przy użyciu - używaj na własne ryzyko.

Przykładowe użycie:

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

Pełny przykład z testami jednostkowymi @ 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;
  }
})();

Implementacja odpowiedzi Seana w 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));

Jestem pewien, że to prawda. Wygląda na to że działa :-)

Te wielokąty wyglądają tak, jeśli zastanawiasz się:

Jest to zaimplementowana funkcja dla OpenLayers 2 . Warunkiem posiadania wielokąta zgodnego z ruchem wskazówek zegara jest area < 0to potwierdzone przez to odniesienie .

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

Jeśli używasz Matlaba, funkcja ispolycwzwraca true, jeśli wierzchołki wielokąta są w kolejności zgodnej z ruchem wskazówek zegara.

Jak wyjaśniono również w tym artykule Wikipedia orientacji Curve , biorąc pod uwagę 3 punkty p, qa rna płaszczyźnie (tj z X i Y), można obliczyć znak następującym wyznacznika

Jeśli wyznacznik jest ujemny (tj. Orient(p, q, r) < 0), Wówczas wielokąt jest zorientowany zgodnie z ruchem wskazówek zegara (CW). Jeśli wyznacznik jest dodatni (tj. Orient(p, q, r) > 0), Wielokąt jest zorientowany przeciwnie do ruchu wskazówek zegara (CCW). Wyznacznik wynosi zero (tzn. Orient(p, q, r) == 0Jeśli punkty p) qi rsą współliniowe .

We wzorze powyżej, dołączana te przed współrzędnych p, q i rdlatego korzysta jednorodne współrzędnych .

Myślę, że aby niektóre punkty były podawane zgodnie z ruchem wskazówek zegara, wszystkie krawędzie muszą być dodatnie, nie tylko suma krawędzi. Jeśli jedna krawędź jest ujemna, co najmniej 3 punkty zostaną podane przeciwnie do ruchu wskazówek zegara.

Moje rozwiązanie C # / LINQ opiera się na poradach dotyczących różnych produktów @charlesbretana poniżej. Możesz określić normalną wartość odniesienia dla uzwojenia. Powinno działać, dopóki krzywa znajduje się głównie w płaszczyźnie określonej przez wektor w górę.

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;

        } 

    }
}

z testem jednostkowym

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



        } 
    }
}

Oto moje rozwiązanie, korzystając z wyjaśnień w innych odpowiedziach:

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

Znacznie prostsza obliczeniowo metoda, jeśli znasz już punkt wewnątrz wielokąta :

1. Wybierz dowolny segment linii z oryginalnego wielokąta, punktów i ich współrzędnych w tej kolejności.

2. Dodaj znany punkt „wewnętrzny” i utwórz trójkąt.

3. Oblicz CW lub CCW zgodnie z sugestią tutaj z tymi trzema punktami.

Po przetestowaniu kilku niewiarygodnych implementacji, algorytmem, który dostarczył zadowalające wyniki dotyczące orientacji CW / CCW po wyjęciu z pudełka, był ten opublikowany przez OP w tym wątku (shoelace_formula_3 ).

Jak zawsze liczba dodatnia reprezentuje orientację CW, podczas gdy liczba ujemna CCW.

Oto szybkie rozwiązanie 3.0 oparte na powyższych odpowiedziach:

    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

Kolejne rozwiązanie tego;

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;
}

Weź wszystkie wierzchołki jako tablicę taką jak ta;

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

Rozwiązanie dla R, aby określić kierunek i odwrócić, jeśli jest zgodny z ruchem wskazówek zegara (uznał to za konieczne dla owin obiektów):

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

Chociaż te odpowiedzi są poprawne, są bardziej matematycznie intensywne niż to konieczne. Załóżmy współrzędne mapy, gdzie najbardziej wysunięty na północ punkt jest najwyższym punktem na mapie. Znajdź punkt na północy, a jeśli 2 punkty wiążą się, to jest to najbardziej na północ, a następnie na wschód (jest to punkt, który lhf używa w swojej odpowiedzi). W twoich punktach

punkt [0] = (5,0)

punkt [1] = (6,4)

punkt [2] = (4,5)

punkt [3] = (1,5)

punkt [4] = (1,0)

Jeśli założymy, że P2 jest najbardziej na północ, a następnie na wschód, to poprzedni lub następny punkt określ zgodnie z ruchem wskazówek zegara, CW lub CCW. Ponieważ najbardziej wysunięty na północ punkt znajduje się na północnej ścianie, jeśli P1 (poprzedni) do P2 przesuwa się na wschód, kierunek to CW. W tym przypadku porusza się na zachód, więc zgodnie z przyjętą odpowiedzią kierunek to CCW. Jeśli poprzedni punkt nie ma ruchu poziomego, wówczas ten sam system dotyczy następnego punktu, P3. Jeśli P3 jest na zachód od P2, to znaczy, że ruch jest w kierunku przeciwnym do ruchu wskazówek zegara. Jeśli ruch P2 do P3 odbywa się na wschód, w tym przypadku na zachód, ruch jest CW. Załóżmy, że nte, P2 w twoich danych, jest najbardziej na północ niż wschodni punkt, a prv to poprzedni punkt, P1 w twoich danych, a nxt to następny punkt, P3 w twoich danych, a [0] jest poziomy lub wschodni / zachód, gdzie zachód jest mniejszy niż wschód, a [1] jest pionowy.

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)

Kod C # do implementacji odpowiedzi lhf :

// 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;
}

Oto prosta implementacja języka Python 3 oparta na tej odpowiedzi (która z kolei oparta jest na rozwiązaniu zaproponowanym w zaakceptowanej odpowiedzi )

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

znajdź środek masy tych punktów.

załóżmy, że są linie od tego punktu do twoich punktów.

znajdź kąt między dwiema liniami dla linii 0 linia1

niż zrób to dla linii 1 i linii 2

...

...

jeśli kąt ten rośnie monotonicznie niż w kierunku przeciwnym do ruchu wskazówek zegara,

w przeciwnym razie monotonicznie zmniejsza się zgodnie z ruchem wskazówek zegara

inaczej (to nie jest monotonne)

nie możesz zdecydować, więc nie jest to mądre