¿Cómo determinar si una lista de puntos de polígono está en el sentido de las agujas del reloj?
https://stackoverflow.com/questions/1165647
-
19-09-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianPregunta
Tener 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).
Solución
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 de ininterroncación como una figura ocho, diciéndole si es así. principalmente agujas del reloj).
Suma sobre los bordes (x2 - x1) (y2 + y1). 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
Otros consejos
los producto cruzado mide el grado de perpendicular de dos vectores. Imagine que cada borde de su polígono es un vector en el plano XY de un espacio XYZ tridimensional (3-D). Luego, el producto cruzado de dos bordes sucesivos es un vector en la dirección Z (dirección Z positiva si el segundo segmento es en sentido horario, menos la 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 se apaga 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)
Así que etiqueta los bordes consecutivamente como
edgeA es el segmento de point0 a point1 y
edgeB Entre point1 a point2
...
edgeE está entre point4 y point0.
Entonces el vértice A (point0) está entre
edgeE De point4 a point0]
edgeA De point0 a `Point1 '
Estos dos bordes son en sí mismos vectores, cuyas coordenadas X e Y se pueden determinar restando las coordenadas de sus puntos de inicio y finalización:
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 adyacentes se calcula utilizando el determinante de la siguiente matriz, que se construye colocando las coordenadas de los dos vectores debajo de los símbolos que representan el eje de las tres coordenadas (i, j, & k). La tercera coordenada valorada (cero) está allí porque el concepto de producto cruzado es una construcción tridimensional, por lo que extendemos 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 k, (o eje z) componente.
La fórmula para calcular la magnitud del k o el componente 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 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 la derecha, y por cuánto. No hay necesidad de tomar el arco. ¡Todo lo 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 agregue 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 tirar otra solución, porque es sencillo y no matemáticamente intensivo, solo usa álgebra básica. Calcule el área firmada del polígono. Si es negativo, los puntos están en sentido horario, si es positivo, son en sentido antihorario. (Esto es muy similar a la solución de Beta).
Calcule el área firmada: a = 1/2 * (x1*y2 - X2*y1 + x2*y3 - X3*y2 + ... + xnorte*y1 - X1*ynorte)
O en el código de pseudo:
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á revisando el pedido, no necesita molestarse en dividir entre 2.
Encuentre el vértice con más pequeño y (y más grande X si hay lazos). Deja que el vértice sea A y el vértice anterior en la lista sea B y el siguiente vértice en la lista sea C. Ahora calcula el señal del producto cruzado de AB y AC.
Referencias:
¿Cómo encuentro la orientación de un polígono simple? enPreguntas frecuentes: comp.graphics.algorithms.
Orientación de la curva en Wikipedia.
Aquí hay una simple implementación de C# del algoritmo basado en esta respuesta.
Supongamos que tenemos un Vector escribir tener X y Y Propiedades del 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 módulo o el operador del resto realizando 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á cero (así que omitalos); Para el resto, el seno del ángulo será dado por el producto cruzado de las normalizaciones a la longitud de la unidad de (punto [N] -Point [0]) y (punto [N-1] -Point [0]).
Si la suma de los valores es positiva, entonces su polígono se dibuja en sentido antihorario.
Para lo que vale, utilicé esta mezcla para calcular el orden de devanado para las aplicaciones API V3 de Google Maps.
El código aprovecha el efecto secundario de las áreas de polígono: un orden de devanado en sentido horario de vértices 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: use 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 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;
}
})();
Una implementación de 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 lo estás preguntando:
Si usa MATLAB, la función ispolycw Devuelve verdadero si los vértices de polígono están en el sentido de las agujas del reloj.
Esta es la función implementada para OpenLayers 2. La condición para tener un polígono en sentido horario es area < 0, confirmado 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);
}
Como también se explica en este artículo de Wikipedia Orientación de la curva, dado 3 puntos p, q y r En 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 los puntos p, q y r son colineal.
En la fórmula anterior, prependemos los que están frente a las coordenadas de p, q y r Porque estamos usando coordenadas homogéneas.
Creo que para que algunos puntos reciban en sentido horario, todos los bordes deben ser positivos no solo la suma de los bordes. Si un borde es negativo que al menos 3 puntos se dan en sentido antihorario.
Mi solución C# / LINQ se basa en el consejo de productos cruzados de @charlesbretana está a continuación. Puede especificar una referencia normal para el devanado. Debería funcionar mientras 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 computacionalmente más simple, Si ya conoces un punto dentro del polígono:
Elija cualquier segmento de línea del polígono original, puntos y sus coordenadas en ese orden.
Agregue un punto "interno" conocido y forme un triángulo.
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, 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;
}
Tome 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 R para determinar la dirección y revertir si en sentido horario (lo encontró necesario para los objetos Owin):
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 más intensas matemáticamente de lo necesario. Suponga que las coordenadas del mapa, donde la mayor cantidad de punto norte es el punto más alto del mapa. Encuentre el punto más norte, y si 2 puntos se unen, es la mayoría del norte que la mayoría del 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 la mayor cantidad de North, East Point, ya sea el anterior o el siguiente punto, determine en sentido horario, CW o CCW. Dado que el punto más 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, lo es, entonces el movimiento es CCW. Si el movimiento P2 a P3 es este, es oeste en este caso, el movimiento es CW. Suponga que NTE, P2 en sus datos, es la mayor cantidad de puntos del Norte y East y el PRV es el punto anterior, P1 en sus datos, y NXT es el siguiente punto, P3 en sus datos, y [0] es horizontal o este// Oeste donde el oeste es menor que el 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 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;
}
Encuentra el centro de masa de estos puntos.
Supongamos que hay líneas desde este punto hasta sus puntos.
Encuentre el ángulo entre dos líneas para Line0 Line1
que hacerlo para Line1 y Line2
...
...
Si este ángulo está aumentando monotónicamente que en sentido antihorario,
de lo contrario, si la disminución monotónica es en sentido horario
más (no es monotonical)
no puedes decidir, así que no es sabio
