¿Cómo determinar si una lista de puntos poligonales está en el sentido de las agujas del reloj?

Teniendo una lista de puntos, ¿cómo encuentro si están en el sentido de las agujas del reloj?

Por ejemplo:

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

diría que es en sentido antihorario (o en sentido antihorario, para algunas personas).

Algunos de los métodos sugeridos fallarán en el caso de un polígono no convexo, como una media luna. Aquí hay uno simple que funcionará con polígonos no convexos (incluso funcionará con un polígono auto intersectante como una figura ocho, que le indica si es principalmente en sentido horario).

Suma sobre los bordes, (x ₂ - x ₁ ) (y ₂ + y ₁ ). Si el resultado es positivo, la curva es en sentido horario, si es negativa, la curva es en sentido antihorario. (El resultado es el doble del área cerrada, con una convención +/-).

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

El producto cruzado mide el grado de perpendicularidad de dos vectores. Imagine que cada borde de su polígono es un vector en el plano xy de un espacio tridimensional (3-D) xyz. Entonces, el producto cruzado de dos aristas sucesivas es un vector en la dirección z, (dirección z positiva si el segundo segmento es en sentido horario, menos dirección z si es en sentido antihorario). La magnitud de este vector es proporcional al seno del ángulo entre los dos bordes originales, por lo que alcanza un máximo cuando son perpendiculares y disminuye gradualmente para desaparecer cuando los bordes son colineales (paralelos).

Entonces, para cada vértice (punto) del polígono, calcule la magnitud del producto cruzado de los dos bordes adyacentes:

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

Por lo tanto, etiquete los bordes consecutivamente como
edgeAes el segmento de point0a point1y
edgeBentre point1a point2
...
edgeEestá entre point4y point0.

Entonces Vertex A ( point0) está entre
edgeE[De point4a point0]
edgeA[De point0a `punto1 '

Estos dos bordes son en sí mismos vectores, cuyas coordenadas x e y se pueden determinar restando las coordenadas de sus puntos inicial y final:

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

Y el producto cruzado de estos dos bordes contiguos se calcula utilizando el determinante de la matriz siguiente, que se construye poniendo las coordenadas de los dos vectores por debajo de los símbolos que representan los tres ejes de coordenadas ( i, j, y k). La tercera coordenada (cero) está ahí porque el concepto de producto cruzado es una construcción tridimensional, por lo que ampliamos estos vectores 2-D a 3-D para aplicar el producto cruzado:

 i    j    k 
-4    0    0
 1    4    0    

Dado que todos los productos cruzados producen un vector perpendicular al plano de dos vectores que se multiplican, el determinante de la matriz anterior solo tiene un kcomponente (o eje z).
La fórmula para calcular la magnitud del kcomponente del eje z es
a1*b2 - a2*b1 = -4* 4 - 0* 1 = -16

La magnitud de este valor ( -16), es una medida del seno del ángulo entre los 2 vectores originales, multiplicado por el producto de las magnitudes de los 2 vectores.
En realidad, otra fórmula para su valor es
A X B (Cross Product) = |A| * |B| * sin(AB).

Entonces, para volver solo a una medida del ángulo, necesita dividir este valor, ( -16), por el producto de las magnitudes de los dos vectores.

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

Entonces la medida del pecado (AB) = -16 / 16.4924=-.97014...

Esta es una medida de si el siguiente segmento después del vértice se ha doblado hacia la izquierda o hacia la derecha, y por cuánto. No hay necesidad de tomar arco-seno. ¡Lo único que nos importará es su magnitud y, por supuesto, su signo (positivo o negativo)!

Haga esto para cada uno de los otros 4 puntos alrededor de la ruta cerrada, y sume los valores de este cálculo en cada vértice.

Si la suma final es positiva, fue en sentido horario, negativo, en sentido antihorario.

Supongo que esta es una pregunta bastante antigua, pero de todos modos voy a lanzar otra solución, porque es sencilla y no es matemáticamente intensiva, solo usa álgebra básica. Calcule el área firmada del polígono. Si es negativo, los puntos están en el sentido de las agujas del reloj, si es positivo, son en sentido contrario. (Esto es muy similar a la solución de Beta).

Calcule el área con signo: A = 1/2 * (x ₁ * y ₂ - x ₂ * y ₁ + x ₂ * y ₃ - x ₃ * y ₂ + ... + x _n * y ₁ - x ₁ * y _n )

O en seudocódigo:

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

Tenga en cuenta que si solo está verificando el pedido, no necesita molestarse en dividir por 2.

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

Encuentra el vértice con la y más pequeña (y la x más grande si hay vínculos). Deje que el vértice sea Ay el vértice anterior en la lista sea By el siguiente vértice en la lista sea C. Ahora calcule el signo del producto cruzado de ABy AC.

Referencias

• ¿Cómo encuentro la orientación de un polígono simple? en Preguntas frecuentes: comp.graphics.algorithms .

• Orientación de curvas en Wikipedia.

Aquí hay una implementación simple de C # del algoritmo basada en esta respuesta .

Supongamos que tenemos un Vectortipo que tiene Xy Ypropiedades de tipo 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;
}

%es el operador de módulo o resto que realiza la operación de módulo que ( según Wikipedia ) encuentra el resto después de la división de un número por otro.

Comience en uno de los vértices y calcule el ángulo subtendido por cada lado.

El primero y el último serán cero (así que omítalos); para el resto, el seno del ángulo estará dado por el producto cruzado de las normalizaciones a la unidad de longitud de (punto [n] -punto [0]) y (punto [n-1] -punto [0]).

Si la suma de los valores es positiva, su polígono se dibuja en sentido antihorario.

Para lo que vale, utilicé este mixin para calcular el orden de liquidación de las aplicaciones de Google Maps API v3.

El código aprovecha el efecto secundario de las áreas de polígono: un orden de vértices en el sentido de las agujas del reloj produce un área positiva, mientras que un orden de devanado en sentido antihorario de los mismos vértices produce la misma área que un valor negativo. El código también utiliza una especie de API privada en la biblioteca de geometría de Google Maps. Me sentí cómodo usándolo, úselo bajo su propio riesgo.

Uso de la muestra:

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

Ejemplo completo con pruebas unitarias @ 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;
  }
})();

Una implementación de la respuesta 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));

Estoy bastante seguro de que esto es correcto. Parece estar funcionando :-)

Esos polígonos se ven así, si te estás preguntando:

Esta es la función implementada para OpenLayers 2 . La condición para tener un polígono en sentido horario es area < 0confirmada por esta referencia .

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 usa Matlab, la función ispolycwdevuelve verdadero si los vértices del polígono están en el orden de las agujas del reloj.

Como también se explica en este artículo de Wikipedia Orientación de curvas , dados 3 puntos p, qy ren el plano (es decir, con coordenadas x e y), puede calcular el signo del siguiente determinante

Si el determinante es negativo (es decir Orient(p, q, r) < 0), entonces el polígono está orientado en sentido horario (CW). Si el determinante es positivo (es decir Orient(p, q, r) > 0), el polígono está orientado en sentido antihorario (CCW). El determinante es cero (es decir Orient(p, q, r) == 0) si puntos p, qy rson colineales .

En la fórmula anterior, anteponemos las que están delante de las coordenadas de p, q y rporque estamos usando coordenadas homogéneas .

Creo que para que algunos puntos se den en el sentido de las agujas del reloj, todos los bordes deben ser positivos, no solo la suma de los bordes. Si un borde es negativo, al menos 3 puntos se dan en sentido antihorario.

Mi solución C # / LINQ se basa en los consejos de productos cruzados de @charlesbretana que se encuentran a continuación. Puede especificar una referencia normal para el devanado. Debería funcionar siempre que la curva esté principalmente en el plano definido por el vector ascendente.

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;

        } 

    }
}

con una prueba unitaria

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



        } 
    }
}

Esta es mi solución usando las explicaciones en las otras respuestas:

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

Un método mucho más simple desde el punto de vista computacional, si ya conoce un punto dentro del polígono :

1. Elija cualquier segmento de línea del polígono original, los puntos y sus coordenadas en ese orden.

2. Agregue un punto "interno" conocido y forme un triángulo.

3. Calcule CW o CCW como se sugiere aquí con esos tres puntos.

Después de probar varias implementaciones poco confiables, el algoritmo que proporcionó resultados satisfactorios con respecto a la orientación CW / CCW fuera de la caja fue el publicado por OP en este hilo (shoelace_formula_3 ).

Como siempre, un número positivo representa una orientación CW, mientras que un número negativo CCW.

Aquí está la solución swift 3.0 basada en las respuestas anteriores:

    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

Otra solución para esto;

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

Toma todos los vértices como una matriz como esta;

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

Solución para que R determine la dirección e invierta si es en el sentido de las agujas del reloj (lo encontró necesario para los objetos):

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

Si bien estas respuestas son correctas, son matemáticamente más intensas de lo necesario. Asuma las coordenadas del mapa, donde el punto más al norte es el punto más alto del mapa. Encuentre el punto más al norte, y si empatan 2 puntos, es el más al norte que el más al este (este es el punto que lhf usa en su respuesta). En tus puntos

punto [0] = (5,0)

punto [1] = (6,4)

punto [2] = (4,5)

punto [3] = (1,5)

punto [4] = (1,0)

Si suponemos que P2 es el punto más al norte, entonces este, ya sea el punto anterior o el siguiente, determina en sentido horario, CW o CCW. Como el punto más al norte está en la cara norte, si el P1 (anterior) a P2 se mueve hacia el este, la dirección es CW. En este caso, se mueve hacia el oeste, por lo que la dirección es CCW como dice la respuesta aceptada. Si el punto anterior no tiene movimiento horizontal, entonces el mismo sistema se aplica al siguiente punto, P3. Si P3 está al oeste de P2, entonces el movimiento es CCW. Si el movimiento P2 a P3 es hacia el este, en este caso es hacia el oeste, el movimiento es CW. Suponga que nte, P2 en sus datos, es el punto más al norte que este y que prv es el punto anterior, P1 en sus datos, y nxt es el siguiente punto, P3 en sus datos, y [0] es horizontal u este / oeste donde oeste es menor que este, y [1] es 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)

Código C # para implementar la respuesta 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;
}

Aquí hay una implementación simple de Python 3 basada en esta respuesta (que, a su vez, se basa en la solución propuesta en la respuesta aceptada )

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

encuentra el centro de masa de estos puntos.

supongamos que hay líneas desde este punto hasta tus puntos.

encuentra el ángulo entre dos líneas para line0 line1

que hacerlo para line1 y line2

...

...

si este ángulo está aumentando monotónicamente de lo contrario,

de lo contrario si disminuye monotónicamente es en sentido horario

de lo contrario (no es monotónico)

no puedes decidir, entonces no es sabio