Comment déterminer si une liste de points de polygone est dans le sens horaire?

Ayant une liste de points, comment puis-je trouver s'ils sont dans le sens horaire?

Par exemple:

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

dirait qu'il est anti-horaire (ou anti-horaire, pour certaines personnes).

Certaines des méthodes suggérées échoueront dans le cas d'un polygone non convexe, comme un croissant. En voici un simple qui fonctionnera avec des polygones non convexes (il fonctionnera même avec un polygone auto-intersecté comme un chiffre huit, vous indiquant s'il est principalement dans le sens horaire).

Somme sur les bords, (x ₂ - x ₁ ) (y ₂ + y ₁ ). Si le résultat est positif, la courbe est dans le sens horaire, s'il est négatif, la courbe est dans le sens antihoraire. (Le résultat est le double de la zone fermée, avec une 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

Le produit croisé mesure le degré de perpendiculaire de deux vecteurs. Imaginez que chaque arête de votre polygone soit un vecteur dans le plan xy d'un espace xyz tridimensionnel (3-D). Alors le produit croisé de deux bords successifs est un vecteur dans la direction z (direction z positive si le deuxième segment est dans le sens horaire, moins direction z si c'est dans le sens antihoraire). La magnitude de ce vecteur est proportionnelle au sinus de l'angle entre les deux bords d'origine, il atteint donc un maximum quand ils sont perpendiculaires, et diminue progressivement pour disparaître lorsque les bords sont colinéaires (parallèles).

Donc, pour chaque sommet (point) du polygone, calculez la magnitude du produit croisé des deux bords adjacents:

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

Étiquetez donc les bords consécutivement, de même que
edgeAle segment de point0vers point1et
edgeBentre point1à point2
...
edgeEest entre point4et point0.

Alors le sommet A ( point0) est entre
edgeE[De point4à point0]
edgeA[De point0à `point1 '

Ces deux arêtes sont elles-mêmes vecteurs, dont les coordonnées x et y peuvent être déterminées en soustrayant les coordonnées de leurs points de début et de fin:

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

Et le produit vectoriel de ces deux bords adjacents est calculée en utilisant le facteur déterminant de la matrice suivante, qui est réalisée en mettant les coordonnées des deux vecteurs ci - dessous les symboles représentant les trois axes de coordonnées ( i, j, et k). La troisième coordonnée (zéro) est là parce que le concept de produit croisé est une construction 3D, et donc nous étendons ces vecteurs 2D en 3D afin d'appliquer le produit croisé:

 i    j    k 
-4    0    0
 1    4    0    

Étant donné que tous les produits croisés produisent un vecteur perpendiculaire au plan de deux vecteurs en cours de multiplication, le déterminant de la matrice ci-dessus n'a qu'une kcomposante, (ou axe z).
La formule pour calculer la magnitude de la kcomposante de l'axe ou z est
a1*b2 - a2*b1 = -4* 4 - 0* 1 = -16

La magnitude de cette valeur ( -16), est une mesure du sinus de l'angle entre les 2 vecteurs d'origine, multipliée par le produit des magnitudes des 2 vecteurs.
En fait, une autre formule pour sa valeur est
A X B (Cross Product) = |A| * |B| * sin(AB).

Donc, pour revenir à une mesure de l'angle, vous devez diviser cette valeur, ( -16), par le produit des grandeurs des deux vecteurs.

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

Donc, la mesure du péché (AB) = -16 / 16.4924=-.97014...

Il s'agit de mesurer si le segment suivant après le sommet s'est courbé à gauche ou à droite, et de combien. Il n'est pas nécessaire de prendre l'arc sinus. Tout ce qui nous importera, c'est son ampleur, et bien sûr son signe (positif ou négatif)!

Faites cela pour chacun des 4 autres points autour du chemin fermé et additionnez les valeurs de ce calcul à chaque sommet.

Si la somme finale est positive, vous êtes allé dans le sens horaire, négatif, anti-horaire.

Je suppose que c'est une assez vieille question, mais je vais quand même jeter une autre solution, car elle est simple et peu mathématique - elle utilise simplement l'algèbre de base. Calculez la zone signée du polygone. S'il est négatif, les points sont dans le sens horaire, s'il est positif, ils sont dans le sens antihoraire. (Ceci est très similaire à la solution Beta.)

Calculez la zone signée: A = 1/2 * (x ₁ * y ₂ - x ₂ * y ₁ + x ₂ * y ₃ - x ₃ * y ₂ + ... + x _n * y ₁ - x ₁ * y _n )

Ou en 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

Notez que si vous ne faites que vérifier la commande, vous n'avez pas à vous soucier de diviser par 2.

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

Trouvez le sommet avec le plus petit y (et le plus grand x s'il y a des liens). Soit le sommet Aet le sommet précédent de la liste Bet le sommet suivant de la liste C. Calculez maintenant le signe du produit croisé de ABet AC.

Références:

• Comment trouver l'orientation d'un polygone simple? dans la foire aux questions: comp.graphics.algorithms .

• Orientation des courbes sur Wikipedia.

Voici une implémentation C # simple de l'algorithme basé sur cette réponse .

Supposons que nous avons un Vectortype ayant Xet des Ypropriétés de 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;
}

%est l'opérateur modulo ou reste effectuant l'opération modulo qui ( selon Wikipedia ) trouve le reste après division d'un nombre par un autre.

Commencez à l'un des sommets et calculez l'angle sous-tendu par chaque côté.

Le premier et le dernier seront nuls (alors sautez-les); pour le reste, le sinus de l'angle sera donné par le produit croisé des normalisations à la longueur unitaire de (point [n] -point [0]) et (point [n-1] -point [0]).

Si la somme des valeurs est positive, votre polygone est dessiné dans le sens inverse des aiguilles d'une montre.

Pour ce que cela vaut, j'ai utilisé ce mixin pour calculer l'ordre d'enroulement pour les applications Google Maps API v3.

Le code exploite l'effet secondaire des zones de polygones: un ordre d'enroulement dans le sens horaire des sommets donne une zone positive, tandis qu'un ordre d'enroulement dans le sens antihoraire des mêmes sommets produit la même zone qu'une valeur négative. Le code utilise également une sorte d'API privée dans la bibliothèque de géométrie de Google Maps. Je me sentais à l'aise de l'utiliser - à vos risques et périls.

Exemple d'utilisation:

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

Exemple complet avec tests unitaires @ 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  [email protected]
 *
 *  @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;
  }
})();

Une implémentation de la réponse de Sean en 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));

Je suis sûr que c'est vrai. Cela semble fonctionner :-)

Ces polygones ressemblent à ceci, si vous vous demandez:

Il s'agit de la fonction implémentée pour OpenLayers 2 . La condition pour avoir un polygone dans le sens horaire est area < 0, cela est confirmé par cette référence .

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

Si vous utilisez Matlab, la fonction ispolycwrenvoie true si les sommets du polygone sont dans le sens des aiguilles d'une montre.

Comme expliqué également dans cet article de Wikipedia Orientation de la courbe , étant donné 3 points p, qet rsur le plan (c'est-à-dire avec les coordonnées x et y), vous pouvez calculer le signe du déterminant suivant

Si le déterminant est négatif (c'est-à-dire Orient(p, q, r) < 0), alors le polygone est orienté dans le sens horaire (CW). Si le déterminant est positif (c'est-à-dire Orient(p, q, r) > 0), le polygone est orienté dans le sens antihoraire (CCW). Le déterminant est zéro (c'est-à-dire Orient(p, q, r) == 0) si points p, qet rsont colinéaires .

Dans la formule ci-dessus, nous ajoutons ceux devant les coordonnées de p, q et rparce que nous utilisons des coordonnées homogènes .

Code C # pour implémenter la réponse de 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;
}

Je pense que pour que certains points soient donnés dans le sens horaire, tous les bords doivent être positifs non seulement la somme des bords. Si un bord est négatif, au moins 3 points sont donnés dans le sens antihoraire.

Ma solution C # / LINQ est basée sur les conseils produits croisés de @charlesbretana ci-dessous. Vous pouvez spécifier une normale de référence pour l'enroulement. Cela devrait fonctionner tant que la courbe est principalement dans le plan défini par le vecteur haut.

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;

        } 

    }
}

avec un test unitaire

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



        } 
    }
}

Voici ma solution en utilisant les explications des autres réponses:

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

Une méthode beaucoup plus simple sur le plan informatique, si vous connaissez déjà un point à l'intérieur du polygone :

1. Choisissez n'importe quel segment de ligne dans le polygone d'origine, les points et leurs coordonnées dans cet ordre.

2. Ajoutez un point «intérieur» connu et formez un triangle.

3. Calculez CW ou CCW comme suggéré ici avec ces trois points.

Après avoir testé plusieurs implémentations peu fiables, l'algorithme qui a fourni des résultats satisfaisants concernant l'orientation CW / CCW hors de la boîte était celui publié par OP dans ce fil ( shoelace_formula_3).

Comme toujours, un nombre positif représente une orientation CW, tandis qu'un nombre négatif CCW.

Voici la solution swift 3.0 basée sur les réponses ci-dessus:

    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

Une autre solution pour cela;

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

Prenez tous les sommets comme un tableau comme celui-ci;

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 pour R pour déterminer la direction et inverser si dans le sens horaire (jugé nécessaire pour les objets propres):

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

Bien que ces réponses soient correctes, elles sont plus mathématiquement intenses que nécessaires. Supposons les coordonnées de la carte, où le point le plus au nord est le point le plus élevé sur la carte. Trouvez le point le plus au nord, et si 2 points sont à égalité, c'est le plus au nord puis le plus à l'est (c'est le point que lhf utilise dans sa réponse). Dans vos points,

point [0] = (5,0)

point [1] = (6,4)

point [2] = (4,5)

point [3] = (1,5)

point [4] = (1,0)

Si nous supposons que P2 est le point le plus au nord puis à l'est, le point précédent ou suivant détermine dans le sens horaire, CW ou CCW. Puisque le point le plus au nord se trouve sur la face nord, si le P1 (précédent) vers P2 se déplace vers l'est, la direction est CW. Dans ce cas, il se déplace vers l'ouest, donc la direction est CCW comme le dit la réponse acceptée. Si le point précédent n'a pas de mouvement horizontal, alors le même système s'applique au point suivant, P3. Si P3 est à l'ouest de P2, c'est le cas, alors le mouvement est CCW. Si le mouvement P2 à P3 est à l'est, c'est à l'ouest dans ce cas, le mouvement est CW. Supposons que nte, P2 dans vos données, est le point le plus au nord puis à l'est et le prv est le point précédent, P1 dans vos données et nxt est le point suivant, P3 dans vos données et [0] est horizontal ou est / ouest où l'ouest est inférieur à l'est et [1] est 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)

Voici une implémentation Python 3 simple basée sur cette réponse (qui, à son tour, est basée sur la solution proposée dans la réponse acceptée )

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

trouver le centre de masse de ces points.

supposons qu'il y ait des lignes de ce point à vos points.

trouver l'angle entre deux lignes pour line0 line1

que de le faire pour la ligne 1 et la ligne 2

...

...

si cet angle augmente de façon monotone que dans le sens antihoraire,

sinon si la diminution est monotone dans le sens horaire

sinon (ce n'est pas monotone)

vous ne pouvez pas décider, donc ce n'est pas sage