如何确定一个多边形点的列表是顺时针顺序?

有了一个点列表,我如何确定它们是否是顺时针顺序的?

例如:

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

会说它是逆时针的(对某些人来说是逆时针的)

177742 次浏览

从其中一个顶点开始,计算每条边对应的角度。

第一个和最后一个将是零(所以跳过它们);对于其余部分,角度的正弦值将由归一化与(点[n]-点[0])和(点[n-1]-点[0])的单位长度的叉乘给出。

如果这些值的和是正的,那么你的多边形是逆时针方向绘制的。

求出这些点的质心。

假设有直线从这个点到你们的点。

求line0 line1的两条直线夹角

而不是直线1和直线2

...

...

如果这个角是单调递增的,而不是逆时针递增的,

如果是单调递减,则是顺时针递减

Else(它不是单调的)

你不能决定,所以这是不明智的

一些建议的方法在非凸多边形(如新月形)的情况下会失败。这里有一个简单的,可以用于非凸多边形(它甚至可以用于自相交的多边形,如数字8,告诉你它是否是主要是顺时针)。

边求和,(x2−x1)(y2 + y1)。如果结果是正的,曲线是顺时针的,如果结果是负的,曲线是逆时针的。(结果是封闭面积的两倍,采用+/-惯例。)

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

叉积测量两个向量的垂直度。假设多边形的每条边都是三维xyz空间x-y平面上的一个向量。那么两条连续边的叉乘就是一个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是从point0point1的段,并且
edgeBpoint1point2之间
< br >… edgeE介于point4point0之间。< / p > 那么顶点A (point0)在
之间 edgeE [From point4 to point0]
edgeA[从point0到' point1'

.

这两条边本身就是向量,它们的x坐标和y坐标可以通过减去它们的起点和终点的坐标来确定:

< p > edgeE = point0 - point4 = (1, 0) - (5, 0) = (-4, 0)和< br > edgeA = point1 - point0 = (6, 4) - (1, 0) = (5, 4)和< / p >

这两条相邻边的叉乘是用下面矩阵的行列式来计算的,这个矩阵是通过将两个向量的坐标放在表示三个坐标轴的符号下面来构造的(ij, &k)。第三个(零)值坐标在那里,因为叉乘概念是一个三维结构,所以我们将这些二维向量扩展到三维,以便应用叉乘:

 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)。< / p >

所以,为了得到角度的测量值,你需要用这个值(-16)除以两个向量大小的乘积。

__abc0 = __abc1 = __abc2

所以sinab = -16 / 16.4924 = -.97014...

这是一个度量顶点后的下一段是否向左或向右弯曲,以及弯曲的程度。不需要取arcsin函数。我们只关心它的大小,当然还有它的符号(正的还是负的)!

对闭合路径周围的其他4个点都这样做,并将每个顶点的计算值相加。

如果最终的和是正的,就顺时针,负的,逆时针。

找出y最小的顶点(如果有平手,则x最大)。设顶点为A,列表中的前一个顶点为B,列表中的下一个顶点为C。现在计算ABAC的叉乘的标志


引用:

我将提出另一个解决方案,因为它很简单,不需要大量的数学运算,它只是使用了基本的代数。计算多边形的带符号面积。如果是负的,点是顺时针的,如果是正的,点是逆时针的。(这与Beta的解决方案非常相似。)

计算有符号面积: = 1/2 * (x1 * y2 - x2 * y1 + x2 * y3. - x3. * y2 +…+ xn*y1 - x1*yn)

或者在伪代码中:

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

这是我使用其他答案中的解释的解决方案:

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

下面是基于@Beta的回答的算法的简单c#实现。

让我们假设我们有一个Vector类型,具有类型为doubleXY属性。

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

%是模运算符或余数运算符,执行模运算,(根据维基百科)在一个数除以另一个数后求余数。


根据@MichelRouzic评论的优化版本:

double sum = 0.0;
Vector v1 = vertices[vertices.Count - 1]; // or vertices[^1] with
// C# 8.0+ and .NET Core
for (int i = 0; i < vertices.Count; i++) {
Vector v2 = vertices[i];
sum += (v2.X - v1.X) * (v2.Y + v1.Y);
v1 = v2;
}
return sum > 0.0;

这不仅节省了模运算%,还节省了数组索引。

如果使用Matlab,如果多边形顶点按顺时针顺序排列,函数ispolycw将返回true。

为了它的价值,我使用这个mixin来计算谷歌Maps API v3应用程序的缠绕顺序。

该代码利用了多边形区域的副作用:顺时针旋转顺序的顶点产生一个正的区域,而逆时针旋转顺序的相同顶点产生一个负的区域。该代码还使用了谷歌Maps几何库中的一种私有API。我觉得使用它很舒服——使用风险自负。

示例用法:

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

正如在这篇维基百科文章曲线方向中所解释的,给定平面上的3个点pqr(即x和y坐标),你可以计算以下行列式的符号

enter image description here

如果行列式为负(即Orient(p, q, r) < 0),则多边形是顺时针方向(CW)。如果行列式是正的(即Orient(p, q, r) > 0),多边形是逆时针方向的(CCW)。如果点pqr共线的,行列式为零(即Orient(p, q, r) == 0)。

在上面的公式中,我们在pqr的坐标前加上1,因为我们使用的是齐次坐标

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

我认为为了使某些点顺时针方向,所有的边都必须是正的而不仅仅是边的和。如果一条边是负的,则逆时针方向给出至少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();






}
}
}

在测试了几个不可靠的实现后,在CW/CCW方向方面提供令人满意结果的算法是由OP在线程(shoelace_formula_3)中发布的算法。

与往常一样,正数表示CW方向,而负数表示CCW方向。

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

我很确定这是对的。这似乎是有效的:-)

这些多边形看起来是这样的,如果你想知道的话:

< img src = " https://i.imgur.com/JvZmmMH.png " alt = " " >

一个计算上更简单的方法,如果你已经知道多边形内的一个点:

  1. 从原始多边形中选择任意线段,按此顺序选择点及其坐标。

  2. 加上一个已知的“内部”点,形成一个三角形。

  3. 计算CW或CCW建议这里与这三个点。

以下是基于上述答案的swift 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确定方向和反向如果顺时针(发现这是必要的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

虽然这些答案是正确的,但它们在数学上的强度比必要的要大。假设地图坐标,其中最北的点是地图上的最高点。找到最北的点,如果两个点相等,它是最北的,然后是最东的(这是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的西边,那么运动是CCW。如果P2到P3的运动方向是向东,那么在这种情况下,运动方向就是向西。假设数据中的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)

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

下面是一个基于这个答案的简单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