有了一个点列表，我如何确定它们是否是顺时针顺序的?
例如:
point[0] = (5,0)
point[1] = (6,4)
point[2] = (4,5)
point[3] = (1,5)
point[4] = (1,0)
会说它是逆时针的(对某些人来说是逆时针的)
从其中一个顶点开始，计算每条边对应的角度。
第一个和最后一个将是零(所以跳过它们);对于其余部分，角度的正弦值将由归一化与(点[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
The cross product measures the degree of perpendicular-ness of two vectors. Imagine that each edge of your polygon is a vector in the x-y plane of a three-dimensional (3-D) xyz space. Then the cross product of two successive edges is a vector in the z-direction, (positive z-direction if the second segment is clockwise, minus z-direction if it's counter-clockwise). The magnitude of this vector is proportional to the sine of the angle between the two original edges, so it reaches a maximum when they are perpendicular, and tapers off to disappear when the edges are collinear (parallel).
因此，对于多边形的每个顶点(点)，计算两条相邻边的叉乘大小:
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的段 点1到点2之间的edgeB .．. edgeE在point4和point0之间。
那么顶点A (point0)在两者之间 edgeE[从点4到点0] 从point0到' point1'
这两条边本身就是向量，它们的x坐标和y坐标可以通过减去它们的起点和终点的坐标来确定:
edgeE = point0 - point4 = (1,0) - (5,0) = (- 4,0) and edgeA = point1 - point0 = (6,4) - (1,0) = (5,4) and
这两个相邻边的外积是用下面矩阵的行列式来计算的，这个矩阵是通过将两个向量的坐标放在表示三个坐标轴的符号(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(叉乘)= |A| * |B| * sin(AB)。
为了得到角度的大小你需要用这个值(-16)除以两个向量大小的乘积。
|A| * |B| = 4 *根号(17)= 16.4924…
所以sin(AB) = -16 / 16.4924 = -.97014…
这是一个度量顶点后的下一段是否向左或向右弯曲，以及弯曲的程度。不需要取arcsin函数。我们只关心它的大小，当然还有它的符号(正的还是负的)!
对闭合路径周围的其他4个点都这样做，并将每个顶点的计算值相加。
如果最终的和是正的，就顺时针，负的，逆时针。
找出y最小的顶点(如果有平手，则x最大)。假设顶点是A，列表中的前一个顶点是B，列表中的下一个顶点是c。现在计算AB和AC的叉乘的符号。
引用:
如何确定一个简单多边形的方向?在 常见问题:计算机。图形。算法。 维基百科的曲线定位。
我将提出另一个解决方案，因为它很简单，不需要大量的数学运算，它只是使用了基本的代数。计算多边形的带符号面积。如果是负的，点是顺时针的，如果是正的，点是逆时针的。(这与Beta的解决方案非常相似。)
计算带符号的面积: A = 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类型，它的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;
}
%是执行模运算的模运算符或余数运算符，该运算符(根据维基百科)在一个数除以另一个数后求余数。
根据@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;
这不仅节省了模运算%，还节省了数组索引。
测试(参见与@WDUK的讨论)
public static bool IsClockwise(IList<(double X, double Y)> vertices)
{
double sum = 0.0;
var v1 = vertices[^1];
for (int i = 0; i < vertices.Count; i++) {
var v2 = vertices[i];
sum += (v2.X - v1.X) * (v2.Y + v1.Y);
Console.WriteLine($"(({v2.X,2}) - ({v1.X,2})) * (({v2.Y,2}) + ({v1.Y,2})) = {(v2.X - v1.X) * (v2.Y + v1.Y)}");
v1 = v2;
}
Console.WriteLine(sum);
return sum > 0.0;
}
public static void Test()
{
Console.WriteLine(IsClockwise(new[] { (-5.0, -5.0), (-5.0, 5.0), (5.0, 5.0), (5.0, -5.0) }));
// infinity Symbol
//Console.WriteLine(IsClockwise(new[] { (-5.0, -5.0), (-5.0, 5.0), (5.0, -5.0), (5.0, 5.0) }));
}
为了它的价值，我使用这个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个点p, q和r(即x和y坐标)，您可以计算以下行列式的符号
如果行列式为负(即定向(p, q, r) < 0)，则多边形是顺时针方向(CW)。如果行列式为正(即定向(p, q, r) > 0)，则多边形是逆时针方向(CCW)。如果点p, q和r共线，行列式为零(即定向(p, q, r) == 0)。
在上面的公式中，由于我们使用的是齐次坐标，我们将1放在p, q和r的坐标前面。
这是OpenLayers 2的实现函数。有一个顺时针多边形的条件是面积< 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);
}
我的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方向。
Sean的答案的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));
我很确定这是对的。这似乎是有效的:-)
这些多边形看起来是这样的，如果你想知道的话:
一个计算上更简单的方法，如果你已经知道多边形内的一个点:
从原始多边形中选择任意线段，按此顺序选择点及其坐标。 加上一个已知的“内部”点，形成一个三角形。 根据以上三点计算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)
If we assume that P2 is the most north then east point either the previous or next point determine clockwise, CW, or CCW. Since the most north point is on the north face, if the P1 (previous) to P2 moves east the direction is CW. In this case, it moves west, so the direction is CCW as the accepted answer says. If the previous point has no horizontal movement, then the same system applies to the next point, P3. If P3 is west of P2, it is, then the movement is CCW. If the P2 to P3 movement is east, it's west in this case, the movement is CW. Assume that nte, P2 in your data, is the most north then east point and the prv is the previous point, P1 in your data, and nxt is the next point, P3 in your data, and [0] is horizontal or east/west where west is less than east, and [1] is 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#代码实现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
Javascript实现的lhf的答案 (再次强调，这只适用于简单的多边形，即不适用于图8)
let polygon = [ {x:5,y:0}, {x:6,y:4}, {x:4,y:5}, {x:1,y:5}, {x:1,y:0} ] document.body.innerHTML += `Polygon ${polygon.map(p=>`(${p.x}, ${p.y})`).join(", ")} is clockwise? ${isPolygonClockwise(polygon)}` let reversePolygon = [] polygon.forEach(point=>reversePolygon.unshift(point)) document.body.innerHTML += `<br/>Polygon ${reversePolygon.map(p=>`(${p.x}, ${p.y})`).join(", ")} is clockwise? ${isPolygonClockwise(reversePolygon)}` function isPolygonClockwise (polygon) { // From http://www.faqs.org/faqs/graphics/algorithms-faq/ "How do I find the orientation of a simple polygon?" // THIS SOMETIMES FAILS if the polygon is a figure 8, or similar shape where it crosses over itself // Take the lowest point (break ties with the right-most). if (polygon.length < 3) { return true // A single point or two points can't be clockwise/counterclockwise } let previousPoint = polygon[0] let lowestPoint = polygon[1] let nextPoint = polygon[2] polygon.forEach((point, index)=>{ if (point.y > lowestPoint.y || (point.y === lowestPoint.y && point.x > lowestPoint.x)) { // larger y values are lower, in svgs // Break ties with furthest right previousPoint = polygon[(index-1) >= (0) ? (index-1) : (polygon.length-1)] lowestPoint = polygon[index] nextPoint = polygon[(index+1) <= (polygon.length-1) ? (index+1) : (0)] } }) // Check the angle between the previous point, that point, and the next point. // If the angle is less than PI radians, the polygon is clockwise let angle = findAngle(previousPoint, lowestPoint, nextPoint) return angle < Math.PI } function findAngle(A,B,C) { var AB = Math.atan2(B.y-A.y, B.x-A.x); var BC = Math.atan2(C.y-B.y, C.x-B.x); if (AB < 0) AB += Math.PI*2 if (BC < 0) BC += Math.PI*2 return BC-AB; }
对于那些不想“重新发明轮子”的人，我认为值得一提的是，这个检查是在一个名为Shapely (github)的漂亮的Python包中实现的(它基于GEOS C/ c++库):
Shapely is a BSD-licensed Python package for manipulation and analysis of planar geometric objects. It is using the widely deployed open-source geometry library GEOS (the engine of PostGIS, and a port of JTS). Shapely wraps GEOS geometries and operations to provide both a feature rich Geometry interface for singular (scalar) geometries and higher-performance NumPy ufuncs for operations using arrays of geometries. Shapely is not primarily focused on data serialization formats or coordinate systems, but can be readily integrated with packages that are.
来源:https://shapely.readthedocs.io/en/stable/
一个给出OP坐标的小例子:
import numpy as np
from shapely.geometry import Polygon
points = np.array([
(5,0),
(6,4),
(4,5),
(1,5),
(1,0)
])
P = Polygon(points)
这是新构造的多边形:
import matplotlib.pyplot as plt
x,y = P.exterior.coords.xy
plt.plot(x,y)
plt.axis('equal')
plt.grid()
plt.show()
你可以直接使用LinearRing的is_ccw属性来检查多边形是CW还是CCW:
type(P.exterior)
>: shapely.geometry.polygon.LinearRing
P.exterior.is_ccw
>: True
如果颠倒:
points = np.flipud(points)
points
>:
array([[1, 0],
[1, 5],
[4, 5],
[6, 4],
[5, 0]])
P1 = Polygon(points)
P1.exterior.is_ccw
>: True
进一步阅读的文档和参考资料:
shaely is_ccw (github): https://github.com/shapely/shapely/blob/eba985c6e0170ecdd90c83592fd0afa7ae793cb8/shapely/predicates.py#L72-L108 Libgeos (github): https://github.com/libgeos/geos GEOS API参考:https://libgeos.org/doxygen/classgeos_1_1algorithm_1_1Orientation.html#a5af93795969b80f97d7997195974d7c8 GEOS实现(github): https://github.com/libgeos/geos/blob/ab0ce6dafdf7f75ec6d234b6c65bb209037dda17/src/algorithm/Orientation.cpp#L43-L133
