按顺时针顺序排序点?

184
给定一个由x,y点组成的数组,如何按照顺时针(以它们的平均中心点为中心)的顺序对该数组中的点进行排序?我的目标是将这些点传递给线创建函数,最终得到的东西看起来相当“稳定”,尽可能凸出,并且没有线段相交。
值得一提的是,我正在使用Lua,但任何伪代码都将不胜感激。
更新:作为参考,这是基于Ciamej的优秀答案的Lua代码(忽略我的“app”前缀):
function appSortPointsClockwise(points)
    local centerPoint = appGetCenterPointOfPoints(points)
    app.pointsCenterPoint = centerPoint
    table.sort(points, appGetIsLess)
    return points
end

function appGetIsLess(a, b)
    local center = app.pointsCenterPoint

    if a.x >= 0 and b.x < 0 then return true
    elseif a.x == 0 and b.x == 0 then return a.y > b.y
    end

    local det = (a.x - center.x) * (b.y - center.y) - (b.x - center.x) * (a.y - center.y)
    if det < 0 then return true
    elseif det > 0 then return false
    end

    local d1 = (a.x - center.x) * (a.x - center.x) + (a.y - center.y) * (a.y - center.y)
    local d2 = (b.x - center.x) * (b.x - center.x) + (b.y - center.y) * (b.y - center.y)
    return d1 > d2
end

function appGetCenterPointOfPoints(points)
    local pointsSum = {x = 0, y = 0}
    for i = 1, #points do pointsSum.x = pointsSum.x + points[i].x; pointsSum.y = pointsSum.y + points[i].y end
    return {x = pointsSum.x / #points, y = pointsSum.y / #points}
end


2
考虑计算通过该点的径向线的角度,然后按角度排序。 - President James K. Polk
如果你不知道的话,Lua有一个内置函数ipairs(tbl),它可以从1到#tbl迭代tbl的索引和值。因此,对于求和计算,你可以这样做,大多数人认为这看起来更清晰:for _, p in ipairs(points) do pointsSum.x = pointsSum.x + p.x; pointsSum.y = pointsSum.y + p.y end - Ponkadoodle
2
@Wallacoloo 这是非常有争议的。此外,在原始的Lua中,ipairs比数字for循环慢得多。 - Alexander Gladysh
我不得不做一些小的更改才能使它适用于我的情况(只是相对于中心比较两个点)。https://gist.github.com/personalnadir/6624172代码中所有与0的比较似乎都假定点分布在原点周围,而不是任意点。我还认为第一个条件会错误地将点排序到中心点下面。谢谢你的代码,它真的很有帮助! - personalnadir
8个回答

227

首先,计算中心点。 然后,使用任何排序算法对点进行排序,但使用特殊的比较例程来确定一个点是否小于另一个点。

您可以通过以下简单的计算来检查一个点(a)相对于中心点是在左侧还是右侧:

det = (a.x - center.x) * (b.y - center.y) - (b.x - center.x) * (a.y - center.y)

如果结果为零,则它们与中心在同一条直线上,如果结果为正或负,则它们在其中一侧,因此一个点将比另一个点先出现。 使用它,您可以构建一个小于关系来比较点并确定它们在排序数组中应该出现的顺序。但是,您必须定义该顺序的开始位置,我的意思是什么角度将是起始角度(例如x轴的正半部分)。
比较函数的代码可能如下所示:
bool less(point a, point b)
{
    if (a.x - center.x >= 0 && b.x - center.x < 0)
        return true;
    if (a.x - center.x < 0 && b.x - center.x >= 0)
        return false;
    if (a.x - center.x == 0 && b.x - center.x == 0) {
        if (a.y - center.y >= 0 || b.y - center.y >= 0)
            return a.y > b.y;
        return b.y > a.y;
    }

    // compute the cross product of vectors (center -> a) x (center -> b)
    int det = (a.x - center.x) * (b.y - center.y) - (b.x - center.x) * (a.y - center.y);
    if (det < 0)
        return true;
    if (det > 0)
        return false;

    // points a and b are on the same line from the center
    // check which point is closer to the center
    int d1 = (a.x - center.x) * (a.x - center.x) + (a.y - center.y) * (a.y - center.y);
    int d2 = (b.x - center.x) * (b.x - center.x) + (b.y - center.y) * (b.y - center.y);
    return d1 > d2;
}

这将按顺时针方向从12点开始排序点。在同一“小时”内的点将按从中心更远的点开始排序。

如果使用整数类型(在Lua中并不存在),则必须确保det、d1和d2变量是能够容纳执行计算结果的类型。

如果您想要实现尽可能凸的外观,则我想您正在寻找Convex Hull。您可以使用Graham Scan计算它。 在此算法中,您还必须按顺时针(或逆时针)从特殊枢轴点开始对点进行排序。然后,您重复简单的循环步骤,每次检查是否向左或向右添加新点到凸壳中,此检查基于与上述比较函数相同的叉积。

编辑:

添加了一个if语句if (a.y - center.y >= 0 || b.y - center.y >=0)以确保具有x = 0和y < 0的点从更远离中心的点开始排序。如果您不关心同一“小时”内的点的顺序,则可以省略此if语句,并始终返回a.y>b.y

通过添加-center.x-center.y来更正第一个if语句。

添加了第二个if语句(a.x - center.x < 0 && b.x - center.x >= 0)。显然,这是一个疏忽。现在可以重新组织if语句,因为某些检查是冗余的。例如,如果第一个if语句中的第一个条件为false,则第二个if语句的第一个条件必须为true。然而,出于简单起见,我决定保留代码不变。编译器很可能会优化代码并产生相同的结果。


33
没有atan(),没有平方根,甚至没有除法。这是计算机图形学思考的一个很好的例子。尽快剔除所有简单的情况,即使在困难的情况下,也要尽可能地计算少量的内容来知道所需的答案。 - RBerteig
2
如果点集是事先已知的,则仅需要O(n * log n)次比较。如果在此期间要添加点,则需要将它们保留在排序集中,例如平衡二叉搜索树。在这种情况下,添加一个新点需要O(log n)次比较,并且对于涉及极坐标的解决方案来说也完全相同。 Translated text: 如果点集是事先已知的,则仅需要O(n * log n)次比较。如果在此期间要添加点,则需要将它们保留在排序集中,例如平衡二叉搜索树。在这种情况下,添加一个新点需要O(log n)次比较,并且对于涉及极坐标的解决方案来说也完全相同。 - ciamej
2
这里是否遗漏了情况:如果(a.x - center.x < 0 && b.x - center.x >= 0) 返回假; - Tom Martin
2
嘿,这个问题很老了,但是:“这将按顺时针顺序从12点开始排序点。”为什么是12点,我怎么才能把它改成6点呢?有人可以告诉我吗? - Ismoh
1
对于任何寻找更多信息的人,这使用向量之间的叉积。 - SirGuy
显示剩余19条评论

30
你所需要的是一个被称为极坐标的系统。在任何编程语言中,从直角坐标系到极坐标系的转换都很容易实现。公式可以在这个部分找到。
转换为极坐标后,只需按角度theta排序即可。

4
这个方法可以起作用,但它会比回答排序问题所需的计算量更多,存在缺陷。实际上,您并不关心实际角度或径向距离,只关心它们的相对顺序。Ciamej的解决方案更好,因为它避免了除法、平方根和三角函数。 - RBerteig
2
我不确定你对“更好”的标准是什么。例如,将所有点相互比较有点浪费计算资源。三角函数不会让成年人感到害怕,对吧? - Iterator
4
三角函数并不可怕。问题在于计算三角函数的成本很高,并且不需要使用它来确定角度的相对顺序。同样,为了将半径排列顺序,你不需要取平方根。从笛卡尔坐标系完全转换到极坐标系将同时执行反正切和平方根。因此,你的答案是正确的,但在计算机图形学或计算几何学的背景下,这可能不是做这件事的最佳方法。 - RBerteig
我们显然都有同样的想法:每个周期都是宝贵的。 :) - Iterator
2
很高兴能够帮忙。感谢您对使用和速度方面的澄清。我的回答侧重于简单易懂,但@ciamej的回答肯定更加计算效率高。我很高兴这个问题被提出 - 我没有意识到在SO上有一个计算几何社区,我学到了一些东西。我可能会向计算几何社区请教其他问题。 :) - Iterator
显示剩余4条评论

26

对于您的问题,一个有趣的替代方案是寻找旅行商问题(TSP)的近似最小值,即连接所有点的最短路线。如果您的点形成凸形,则应该是正确的解决方案,否则它看起来仍然不错(“固体”形状可以定义为周长/面积比较低,这就是我们在优化的内容)。

您可以使用任何TSP优化器的实现,我相信您在所选择的编程语言中可以找到许多这样的实现。


2
哎呀,"有趣"真是个轻描淡写的说法。 :) - Iterator
1
@迭代器:我对我的想法感到非常满意,但因为它被踩了我感到相当失望 :-/ 你认为这个想法是有效的吗? - static_rtti
我没有给你的建议点踩,但请考虑一下它的计算复杂度。 - Iterator
2
我建议使用其中一种快速近似算法,而不是 NP 完全的原始算法。 - static_rtti
8
感谢提供额外的观点!如果未来有人在寻找备选方案时偶然看到这个帖子,那么具有多个有效但非常不同的答案可能会帮助他们进行头脑风暴。请注意,本翻译力求让内容更加通俗易懂,但并未改变原意。 - Philipp Lenssen
2
请注意,我的方法可能较慢,但在复杂情况下更正确:例如,想象一下“8”的点的情况。在这种情况下,极坐标无法帮助您,而且您得到的结果将严重取决于您选择的中心。TSP解决方案不依赖于任何“启发式”参数。 - static_rtti

3
另一版本(如果a在逆时针方向上出现在b之前,则返回true):
    bool lessCcw(const Vector2D &center, const Vector2D &a, const Vector2D &b) const
    {
        // Computes the quadrant for a and b (0-3):
        //     ^
        //   1 | 0
        //  ---+-->
        //   2 | 3

        const int dax = ((a.x() - center.x()) > 0) ? 1 : 0;
        const int day = ((a.y() - center.y()) > 0) ? 1 : 0;
        const int qa = (1 - dax) + (1 - day) + ((dax & (1 - day)) << 1);

        /* The previous computes the following:

           const int qa =
           (  (a.x() > center.x())
            ? ((a.y() > center.y())
                ? 0 : 3)
            : ((a.y() > center.y())
                ? 1 : 2)); */

        const int dbx = ((b.x() - center.x()) > 0) ? 1 : 0;
        const int dby = ((b.y() - center.y()) > 0) ? 1 : 0;
        const int qb = (1 - dbx) + (1 - dby) + ((dbx & (1 - dby)) << 1);

        if (qa == qb) {
            return (b.x() - center.x()) * (a.y() - center.y()) < (b.y() - center.y()) * (a.x() - center.x());
        } else {
            return qa < qb;
       } 
    }

这样会更快,因为编译器(在Visual C++ 2015上测试)不会生成跳转来计算dax、day、dbx、dby。以下是编译器的输出汇编代码:
; 28   :    const int dax = ((a.x() - center.x()) > 0) ? 1 : 0;

    vmovss  xmm2, DWORD PTR [ecx]
    vmovss  xmm0, DWORD PTR [edx]

; 29   :    const int day = ((a.y() - center.y()) > 0) ? 1 : 0;

    vmovss  xmm1, DWORD PTR [ecx+4]
    vsubss  xmm4, xmm0, xmm2
    vmovss  xmm0, DWORD PTR [edx+4]
    push    ebx
    xor ebx, ebx
    vxorps  xmm3, xmm3, xmm3
    vcomiss xmm4, xmm3
    vsubss  xmm5, xmm0, xmm1
    seta    bl
    xor ecx, ecx
    vcomiss xmm5, xmm3
    push    esi
    seta    cl

; 30   :    const int qa = (1 - dax) + (1 - day) + ((dax & (1 - day)) << 1);

    mov esi, 2
    push    edi
    mov edi, esi

; 31   : 
; 32   :    /* The previous computes the following:
; 33   : 
; 34   :    const int qa =
; 35   :        (   (a.x() > center.x())
; 36   :         ? ((a.y() > center.y()) ? 0 : 3)
; 37   :         : ((a.y() > center.y()) ? 1 : 2));
; 38   :    */
; 39   : 
; 40   :    const int dbx = ((b.x() - center.x()) > 0) ? 1 : 0;

    xor edx, edx
    lea eax, DWORD PTR [ecx+ecx]
    sub edi, eax
    lea eax, DWORD PTR [ebx+ebx]
    and edi, eax
    mov eax, DWORD PTR _b$[esp+8]
    sub edi, ecx
    sub edi, ebx
    add edi, esi
    vmovss  xmm0, DWORD PTR [eax]
    vsubss  xmm2, xmm0, xmm2

; 41   :    const int dby = ((b.y() - center.y()) > 0) ? 1 : 0;

    vmovss  xmm0, DWORD PTR [eax+4]
    vcomiss xmm2, xmm3
    vsubss  xmm0, xmm0, xmm1
    seta    dl
    xor ecx, ecx
    vcomiss xmm0, xmm3
    seta    cl

; 42   :    const int qb = (1 - dbx) + (1 - dby) + ((dbx & (1 - dby)) << 1);

    lea eax, DWORD PTR [ecx+ecx]
    sub esi, eax
    lea eax, DWORD PTR [edx+edx]
    and esi, eax
    sub esi, ecx
    sub esi, edx
    add esi, 2

; 43   : 
; 44   :    if (qa == qb) {

    cmp edi, esi
    jne SHORT $LN37@lessCcw

; 45   :        return (b.x() - center.x()) * (a.y() - center.y()) < (b.y() - center.y()) * (a.x() - center.x());

    vmulss  xmm1, xmm2, xmm5
    vmulss  xmm0, xmm0, xmm4
    xor eax, eax
    pop edi
    vcomiss xmm0, xmm1
    pop esi
    seta    al
    pop ebx

; 46   :    } else {
; 47   :        return qa < qb;
; 48   :    }
; 49   : }

    ret 0
$LN37@lessCcw:
    pop edi
    pop esi
    setl    al
    pop ebx
    ret 0
?lessCcw@@YA_NABVVector2D@@00@Z ENDP            ; lessCcw

享受。


1
switch 中的两个返回语句在数学上是等价的。使用 switch 的原因是什么? - unagi
对于任何寻找 switch(由 @unagi 提到)的人:它已经被编辑删除了。 - Pirulax

2

我知道这篇文章有一个很好的被接受的答案,但是我觉得我还可以贡献一些有用的东西。到目前为止,所有的答案基本上都使用比较函数来比较两个点并确定它们的顺序,但如果你只想使用一个点和一个关键函数怎么办?

不仅如此,而且结果代码也非常紧凑。以下是使用Python内置的sorted函数的完整解决方案:

# Create some random points
num = 7
points = np.random.random((num, 2))
# Compute their center
center = np.mean(points, axis=0)

# Make arctan2 function that returns a value from [0, 2 pi) instead of [-pi, pi)
arctan2 = lambda s, c: angle if (angle := np.arctan2(s, c)) >= 0 else 2 * np.pi + angle

# Define the key function
def clockwise_around_center(point):
    diff = point - center
    rcos = np.dot(diff, center)
    rsin = np.cross(diff, center)
    return arctan2(rsin, rcos)

# Sort our points using the key function
sorted_points = sorted(points, key=clockwise_around_center)

如果这些点在嵌入三维空间的二维平面上,这个答案也适用于三维。我们只需要通过将其与平面的法向量进行点积来修改rsin的计算。例如:

rsin = np.dot([0,0,1], np.cross(diff, center))

如果飞机的法向量是e_z,那么这个代码的优点在于它仅使用一个关键函数逐个处理每个点。如果你在系数级别上计算rsin,你会发现它与接受答案中所称的det完全相同,只是我将其计算在point - centercenter之间,而不是point1 - centerpoint2 - center之间。但是,这个量的几何意义是半径乘以正弦值,因此我将此变量称为rsin。对于点积也是如此,它是半径乘以角度余弦值,因此被称为rcos
有人可能会认为这种解决方案使用了arctan2,因此不够简洁。然而,我个人认为使用关键函数的清晰性比需要调用三角函数更重要。请注意,我更喜欢让arctan2返回[0, 2 pi)中的值,因为当point恰好等于center时,我们得到角度0,因此它将是我们排序列表中的第一个点。这是一个可选项。
为了理解这个代码为什么有效,关键的洞察力在于我们所有的点都是相对于原点定义的箭头,包括center点本身。因此,如果我们计算point - center,这相当于将从center指向point的箭头放置在原点处。因此,我们可以根据它与指向center的箭头所成的角度对箭头point - center进行排序。

1

使用numpy:

import matplotlib.pyplot as plt
import numpy as np

# List of coords
coords = np.array([7,7, 5, 0, 0, 0, 5, 10, 10, 0, 0, 5, 10, 5, 0, 10, 10, 10]).reshape(-1, 2)
centroid = np.mean(coords, axis=0)
sorted_coords = coords[np.argsort(np.arctan2(coords[:, 1] - centroid[1], coords[:, 0] - centroid[0])), :]

plt.scatter(coords[:,0],coords[:,1])
plt.plot(coords[:,0],coords[:,1])
plt.plot(sorted_coords[:,0],sorted_coords[:,1])
plt.show()

enter image description here


1
  • vector3 a = new vector3(1, 0, 0)..............相对于 X 轴
  • vector3 b = 任意点 - 中心点;
- y = |a * b|   ,   x =  a . b

- Atan2(y , x)...............................gives angle between -PI  to  + PI  in radians
- (Input % 360  +  360) % 360................to convert it from  0 to 2PI in radians
- sort by adding_points to list_of_polygon_verts by angle  we got 0  to 360

最终你已经将逆时针排序的顶点解决了

list.Reverse()..................Clockwise_order


0

这里有一种按顺时针排序矩形顶点的方法。我修改了pyimagesearch提供的原始解决方案,并且摆脱了scipy的依赖。

import numpy as np

def pointwise_distance(pts1, pts2):
    """Calculates the distance between pairs of points

    Args:
        pts1 (np.ndarray): array of form [[x1, y1], [x2, y2], ...]
        pts2 (np.ndarray): array of form [[x1, y1], [x2, y2], ...]

    Returns:
        np.array: distances between corresponding points
    """
    dist = np.sqrt(np.sum((pts1 - pts2)**2, axis=1))
    return dist

def order_points(pts):
    """Orders points in form [top left, top right, bottom right, bottom left].
    Source: https://www.pyimagesearch.com/2016/03/21/ordering-coordinates-clockwise-with-python-and-opencv/

    Args:
        pts (np.ndarray): list of points of form [[x1, y1], [x2, y2], [x3, y3], [x4, y4]]

    Returns:
        [type]: [description]
    """
    # sort the points based on their x-coordinates
    x_sorted = pts[np.argsort(pts[:, 0]), :]

    # grab the left-most and right-most points from the sorted
    # x-roodinate points
    left_most = x_sorted[:2, :]
    right_most = x_sorted[2:, :]

    # now, sort the left-most coordinates according to their
    # y-coordinates so we can grab the top-left and bottom-left
    # points, respectively
    left_most = left_most[np.argsort(left_most[:, 1]), :]
    tl, bl = left_most

    # now that we have the top-left coordinate, use it as an
    # anchor to calculate the Euclidean distance between the
    # top-left and right-most points; by the Pythagorean
    # theorem, the point with the largest distance will be
    # our bottom-right point. Note: this is a valid assumption because
    # we are dealing with rectangles only.
    # We need to use this instead of just using min/max to handle the case where
    # there are points that have the same x or y value.
    D = pointwise_distance(np.vstack([tl, tl]), right_most)
    
    br, tr = right_most[np.argsort(D)[::-1], :]

    # return the coordinates in top-left, top-right,
    # bottom-right, and bottom-left order
    return np.array([tl, tr, br, bl], dtype="float32")


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接