Как определить, находится ли список точек многоугольника по часовой стрелке?

Имея список точек, как мне найти, если они расположены по часовой стрелке?

Например:

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

сказал бы, что это против часовой стрелки (или против часовой стрелки, для некоторых людей).

Некоторые из предложенных методов потерпят неудачу в случае невыпуклого многоугольника, такого как полумесяц. Вот простой, который будет работать с невыпуклыми многоугольниками (он даже будет работать с самопересекающимся многоугольником, подобным восьмерке, сообщая, будет ли он в основном по часовой стрелке).

Суммируйте по краям, (x ₂ - x ₁ ) (y ₂ + y ₁ ). Если результат положительный, кривая идет по часовой стрелке, если она отрицательная, кривая против часовой стрелки. (Результат - удвоенная закрытая область с условием +/-.)

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

Крест продукт измеряет степень перпендикулярной-ность двух векторов. Представьте, что каждое ребро вашего многоугольника является вектором в плоскости xy трехмерного (3-D) пространства xyz. Тогда перекрестным произведением двух последовательных ребер является вектор в направлении z (положительное направление z, если второй сегмент направлен по часовой стрелке, минус направление z, если против часовой стрелки). Величина этого вектора пропорциональна синусу угла между двумя исходными краями, поэтому он достигает максимума, когда они перпендикулярны, и сужается, чтобы исчезнуть, когда края коллинеарны (параллельны).

Итак, для каждой вершины (точки) многоугольника рассчитайте величину двойного произведения двух соседних ребер:

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

Так обозначьте края последовательно как
edgeAсегмент от point0к point1и
edgeBмежду point1к point2
...
edgeEмежду point4и point0.

Тогда вершина A ( point0) находится между
edgeE[From point4to point0]
edgeA[From point0to `point1 '

Эти два ребра сами являются векторами, координаты x и y которых можно определить, вычитая координаты их начальной и конечной точек:

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

И крест произведение этих двух смежных ребер вычисляется с использованием определитель следующей матрицы, которая строится путем ввода координат двух векторов ниже символов , представляющих три оси координат ( i, j, & k). Третья (нулевая) -значная координата существует потому, что концепция кросс-произведения является трехмерной конструкцией, и поэтому мы расширяем эти 2-D векторы в 3-D, чтобы применить кросс-произведение:

 i    j    k 
-4    0    0
 1    4    0    

Учитывая, что все перекрестные произведения дают вектор, перпендикулярный к плоскости умножаемых двух векторов, определитель матрицы выше имеет только kкомпонент (или ось Z).
Формула для вычисления величины компонента kили оси z имеет вид
a1*b2 - a2*b1 = -4* 4 - 0* 1 = -16

Величина этого значения ( -16) является мерой синуса угла между двумя исходными векторами, умноженного на произведение величин двух векторов.
На самом деле, другая формула для его стоимости
A X B (Cross Product) = |A| * |B| * sin(AB).

Итак, чтобы вернуться к показателю угла, вам нужно разделить это значение ( -16) на произведение двух векторов.

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

Таким образом, мера греха (AB) = -16 / 16.4924=-.97014...

Это мера того, изогнулся ли следующий сегмент после вершины влево или вправо, и насколько. Нет необходимости принимать арк-синус. Все, о чем мы будем заботиться, это его величина и, конечно, его знак (положительный или отрицательный)!

Сделайте это для каждой из 4 других точек по замкнутому пути и сложите значения из этого расчета в каждой вершине.

Если итоговая сумма положительна, вы пошли по часовой стрелке, отрицательно, против часовой стрелки.

Я думаю, это довольно старый вопрос, но я все равно собираюсь выкинуть другое решение, потому что оно простое и не математически интенсивное - оно просто использует базовую алгебру. Вычислите подписанную площадь многоугольника. Если оно отрицательное, точки расположены по часовой стрелке, если оно положительное, они против часовой стрелки. (Это очень похоже на решение Беты.)

Вычислите область со знаком: A = 1/2 * (x ₁ * y ₂ - x ₂ * y ₁ + x ₂ * y ₃ - x ₃ * y ₂ + ... + x _n * y ₁ - x ₁ * y _н )

Или в псевдокоде:

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

Обратите внимание, что если вы проверяете только порядок, вам не нужно делить на 2.

Источники: http://mathworld.wolfram.com/PolygonArea.html

Найдите вершину с наименьшим y (и наибольшим x, если есть связи). Пусть вершина будет, Aа предыдущая вершина в списке будет, Bа следующая вершина в списке будет C. Теперь вычислите знак перекрестного произведения ABи AC.

Ссылки:

• Как мне найти ориентацию простого многоугольника? в часто задаваемых вопросах: comp.graphics.algorithms .

• Кривая ориентации в Википедии.

Вот простая реализация алгоритма на C #, основанная на этом ответе .

Давайте предположим , что у нас есть Vectorтип , имеющий Xи Yсвойства типа 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;
}

%является оператором по модулю или остатку, выполняющим операцию по модулю, которая ( согласно Википедии ) находит остаток после деления одного числа на другое.

Начните с одной из вершин и вычислите угол, образованный каждой стороной.

Первый и последний будут равны нулю (так что пропустите); в остальном синус угла будет задан как произведение произведенных нормализаций на единицу длины (точка [n] -точка [0]) и (точка [n-1] -точка [0]).

Если сумма значений положительна, то ваш многоугольник нарисован против часовой стрелки.

Для чего бы это ни стоило, я использовал этот миксин для расчета порядка намотки для приложений Google Maps API v3.

В коде используется побочный эффект областей многоугольника: порядок вершин по часовой стрелке дает положительную область, а порядок вращения тех же вершин против часовой стрелки создает ту же площадь, что и отрицательное значение. Код также использует своего рода частный API в библиотеке геометрии Карт Google. Я чувствовал себя комфортно, используя его - используйте на свой страх и риск.

Пример использования:

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

Полный пример с юнит-тестами @ 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;
  }
})();

Реализация ответа Шона в 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));

Уверен, что это правильно. Кажется, это работает :-)

Эти полигоны выглядят так, если вам интересно:

Это реализованная функция для OpenLayers 2 . Условие наличия многоугольника по часовой стрелке area < 0подтверждается этой ссылкой .

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

Если вы используете Matlab, функция ispolycwвозвращает true, если вершины многоугольника расположены по часовой стрелке.

Кроме того, как описано в этой статье Википедии ориентация кривой , давалось 3 очка p, qи rна самолете (т.е. с координатами х и у), можно вычислить знак следующего определителя

Если определитель отрицателен (то есть Orient(p, q, r) < 0), то многоугольник ориентирован по часовой стрелке (CW). Если определитель является положительным (то есть Orient(p, q, r) > 0), многоугольник ориентирован против часовой стрелки (против часовой стрелки). Определитель равен нулю (т.е. Orient(p, q, r) == 0) , если точек p, qи rявляются коллинеарны .

В приведенной выше формуле, мы предварять те , перед координатами p, q и rпотому , что мы используем однородные координаты .

Код C # для реализации ответа 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;
}

Я думаю, что для того, чтобы некоторые точки давались по часовой стрелке, все ребра должны быть положительными, а не только сумма ребер. Если одно ребро отрицательно, то по крайней мере 3 очка даются против часовой стрелки.

Мое решение C # / LINQ основано на совете по продуктам @charlesbretana ниже. Вы можете указать эталонную норму для обмотки. Он должен работать до тех пор, пока кривая в основном находится в плоскости, определенной вектором вверх.

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;

        } 

    }
}

с модульным тестом

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



        } 
    }
}

Это мое решение, используя объяснения в других ответах:

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

Гораздо более простой в вычислительном отношении метод, если вы уже знаете точку внутри многоугольника :

1. Выберите любой отрезок из исходного многоугольника, точки и их координаты в указанном порядке.

2. Добавьте известную «внутреннюю» точку и сформируйте треугольник.

3. Рассчитайте CW или CCW, как предложено здесь, с этими тремя точками.

После тестирования нескольких ненадежных реализаций алгоритм, который дал удовлетворительные результаты в отношении ориентации CW / CCW, был опубликован OP в этой теме ( shoelace_formula_3).

Как всегда, положительное число представляет ориентацию CW, а отрицательное число - CCW.

Вот быстрое решение 3.0, основанное на ответах выше:

    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

Другое решение для этого;

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

Возьмите все вершины в виде такого массива;

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

Решение для R, чтобы определить направление и повернуть, если по часовой стрелке (нашел это необходимым для объектов):

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

Хотя эти ответы верны, они математически более интенсивны, чем необходимо. Предположим, координаты карты, где самая северная точка является самой высокой точкой на карте. Найдите самую северную точку, и если 2 точки связаны, это самый северный, а затем самый восточный (это точка, которую lhf использует в своем ответе). В ваших очках,

точка [0] = (5,0)

точка [1] = (6,4)

точка [2] = (4,5)

точка [3] = (1,5)

точка [4] = (1,0)

Если мы предположим, что P2 - самая северная точка, то восточная точка предыдущей или следующей точки определяет по часовой стрелке, CW или CCW. Так как самая северная точка находится на северной стороне, если P1 (предыдущий) к P2 перемещается на восток, направление - CW. В этом случае он движется на запад, поэтому в направлении принятого ответа указано CCW. Если предыдущая точка не имеет горизонтального перемещения, то та же система применяется к следующей точке, P3. Если P3 к западу от P2, это так, тогда движение - это против часовой стрелки. Если движение P2 к P3 - восток, в этом случае это запад, движение CW. Предположим, что nte, P2 в ваших данных, является самой северной, чем восточной точкой, а prv - предыдущей точкой, P1 в ваших данных, а nxt - следующей точкой, P3 в ваших данных, а [0] - горизонтальная или восточная / запад, где запад меньше, чем восток, а [1] вертикальный.

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)

Вот простая реализация Python 3, основанная на этом ответе (которая, в свою очередь, основана на решении, предложенном в принятом ответе )

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

найти центр масс этих точек.

Предположим, есть линии от этой точки до ваших точек.

найти угол между двумя линиями для line0 line1

чем сделать это для line1 и line2

...

...

если этот угол монотонно увеличивается, чем против часовой стрелки,

иначе, если монотонно уменьшается, то по часовой стрелке

иначе (это не монотонно)

Вы не можете решить, так что это не мудро