Sort points in clockwise order?

Question

Given an array of x,y points, how do I sort the points of this array in clockwise order (around their overall average center point)? My goal is to pass the points to a line-creation function to end up with something looking rather "solid", as convex as possible with no lines intersecting.

For what it's worth, I'm using Lua, but any pseudocode would be appreciated.

Update: For reference, this is the Lua code based on Ciamej's excellent answer (ignore my "app" prefix):

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

Answer 1

First, compute the center point. Then sort the points using whatever sorting algorithm you like, but use special comparison routine to determine whether one point is less than the other.

You can check whether one point (a) is to the left or to the right of the other (b) in relation to the center by this simple calculation:

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

if the result is zero, then they are on the same line from the center, if it's positive or negative, then it is on one side or the other, so one point will precede the other. Using it you can construct a less-than relation to compare points and determine the order in which they should appear in the sorted array. But you have to define where is the beginning of that order, I mean what angle will be the starting one (e.g. the positive half of x-axis).

The code for the comparison function can look like this:

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

This will order the points clockwise starting from the 12 o'clock. Points on the same "hour" will be ordered starting from the ones that are further from the center.

If using integer types (which are not really present in Lua) you'd have to assure that det, d1 and d2 variables are of a type that will be able to hold the result of performed calculations.

If you want to achieve something looking solid, as convex as possible, then I guess you're looking for a Convex Hull. You can compute it using the Graham Scan. In this algorithm, you also have to sort the points clockwise (or counter-clockwise) starting from a special pivot point. Then you repeat simple loop steps each time checking if you turn left or right adding new points to the convex hull, this check is based on a cross product just like in the above comparison function.

Edit:

Added one more if statement if (a.y - center.y >= 0 || b.y - center.y >=0) to make sure that points that have x=0 and negative y are sorted starting from the ones that are further from the center. If you don't care about the order of points on the same 'hour' you can omit this if statement and always return a.y > b.y.

Corrected the first if statements with adding -center.x and -center.y.

Added the second if statement (a.x - center.x < 0 && b.x - center.x >= 0). It was an obvious oversight that it was missing. The if statements could be reorganized now because some checks are redundant. For example, if the first condition in the first if statement is false, then the first condition of the second if must be true. I decided, however, to leave the code as it is for the sake of simplicity. It's quite possible that the compiler will optimize the code and produce the same result anyway.

Answer 2

What you're asking for is a system known as polar coordinates. Conversion from Cartesian to polar coordinates is easily done in any language. The formulas can be found in this section.

After converting to polar coordinates, just sort by the angle, theta.

Answer 3

An interesting alternative approach to your problem would be to find the approximate minimum to the Traveling Salesman Problem (TSP), ie. the shortest route linking all your points. If your points form a convex shape, it should be the right solution, otherwise, it should still look good (a "solid" shape can be defined as one that has a low perimeter/area ratio, which is what we are optimizing here).

You can use any implementation of an optimizer for the TSP, of which I am pretty sure you can find a ton in your language of choice.

Answer 4

Another version (return true if a comes before b in counterclockwise direction):

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

This is faster, because the compiler (tested on Visual C++ 2015) doesn't generate jump to compute dax, day, dbx, dby. Here the output assembly from the compiler:

; 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

Enjoy.

Answer 5

I know this is somewhat of an old post with an excellent accepted answer, but I feel like I can still contribute something useful. All the answers so far essentially use a comparison function to compare two points and determine their order, but what if you want to use only one point at a time and a key function?

Not only is this possible, but the resulting code is also extremely compact. Here is the complete solution using Python's built-in sorted function:

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

This answer would also work in 3D, if the points are on a 2D plane embedded in 3D. We would only have to modify the calculation of rsin by dotting it with the normal vector of the plane. E.g.

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

if that plane has e_z as its normal vector.

The advantage of this code is that it works on only one point at the time using a key function. The quantity rsin, if you work it out on a coefficient level, is exactly the same as what is called det in the accepter answer, except that I compute it between point - center and center, not between point1 - center and point2 - center. But the geometrical meaning of this quantity is the radius times the sin of the angle, hence I call this variable rsin. Similarly for the dot product, which is the radius times the cosine of the angle and hence called rcos.

One could argue that this solution uses arctan2, and is therefore less clean. However, I personally think that the clearity of using a key function outweighs the need for one call to a trig function. Note that I prefer to have arctan2 return a value from [0, 2 pi), because then we get the angle 0 when point happens to be identical to center, and thus it will be the first point in our sorted list. This is an optional choice.

In order to understand why this code works, the crucial insight is that all our points are defined as arrows with respect to the origin, including the center point itself. So if we calculate point - center, this is equivalent to placing the arrow from the tip of center to the tip of point, at the origin. Hence we can sort the arrow point - center with respect to the angle it makes with the arrow pointing to center.

Answer 6

• vector3 a = new vector3(1 , 0 , 0)..............w.r.t X_axis
• vector3 b = any_point - Center;

- 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

Finally you get Anticlockwize sorted verts

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

Answer 7

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

Answer 8

Here's a way to sort the vertices of a rectangle in clock-wise order. I modified the original solution provided by pyimagesearch and got rid of the scipy dependency.

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