Sort 4 3D coordinates in a winding order in any given direction

Question

I need to sort a selection of 3D coordinates in a winding order as seen in the image below. The bottom-right vertex should be the first element of the array and the bottom-left vertex should be the last element of the array. This needs to work given any direction that the camera is facing the points and at any orientation of those points. Since "top-left","bottom-right", etc is relative, I assume I can use the camera as a reference point? We can also assume all 4 points will be coplanar.

I am using the Blender API (writing a Blender plugin) and have access to the camera's view matrix if that is even necessary. Mathematically speaking is this even possible if so how? Maybe I am overcomplicating things?

Since the Blender API is in Python I tagged this as Python, but I am fine with pseudo-code or no code at all. I'm mainly concerned with how to approach this mathematically as I have no idea where to start.

Answer 1

Since you assume the four points are coplanar, all you need to do is find the centroid, calculate the vector from the centroid to each point, and sort the points by the angle of the vector.

import numpy as np

def sort_points(pts):
    centroid = np.sum(pts, axis=0) / pts.shape[0]
    vector_from_centroid = pts - centroid
    vector_angle = np.arctan2(vector_from_centroid[:, 1], vector_from_centroid[:, 0]) 
    sort_order = np.argsort(vector_angle) # Find the indices that give a sorted vector_angle array

    # Apply sort_order to original pts array. 
    # Also returning centroid and angles so I can plot it for illustration. 
    return (pts[sort_order, :], centroid, vector_angle[sort_order])

This function calculates the angle assuming that the points are two-dimensional, but if you have coplanar points then it should be easy enough to find the coordinates in the common plane and eliminate the third coordinate.

Let's write a quick plot function to plot our points:

from matplotlib import pyplot as plt

def plot_points(pts, centroid=None, angles=None, fignum=None):
    fig = plt.figure(fignum)
    plt.plot(pts[:, 0], pts[:, 1], 'or')
    if centroid is not None:
        plt.plot(centroid[0], centroid[1], 'ok')
        
    for i in range(pts.shape[0]):
        lstr = f"pt{i}"
        if angles is not None:
            lstr += f" ang: {angles[i]:.3f}"
        plt.text(pts[i, 0], pts[i, 1], lstr)
    
    return fig

And now let's test this:

With random points:

pts = np.random.random((4, 2))
spts, centroid, angles = sort_points(pts)
plot_points(spts, centroid, angles)

With points in a rectangle:

pts = np.array([[0, 0],  # pt0
                [10, 5], # pt2
                [10, 0], # pt1
                [0, 5]]) # pt3
spts, centroid, angles = sort_points(pts)
plot_points(spts, centroid, angles)

---

It's easy enough to find the normal vector of the plane containing our points, it's simply the (normalized) cross product of the vectors joining two pairs of points:

plane_normal = np.cross(pts[1, :] - pts[0, :], pts[2, :] - pts[0, :])
plane_normal = plane_normal / np.linalg.norm(plane_normal)

Now, to find the projections of all points in this plane, we need to know the "origin" and basis of the new coordinate system in this plane. Let's assume that the first point is the origin, the x axis joins the first point to the second, and since we know the z axis (plane normal) and x axis, we can calculate the y axis.

new_origin = pts[0, :]
new_x = pts[1, :] - pts[0, :]
new_x = new_x / np.linalg.norm(new_x)

new_y = np.cross(plane_normal, new_x)

Now, the projections of the points onto the new plane are given by this answer:

proj_x = np.dot(pts - new_origin, new_x)
proj_y = np.dot(pts - new_origin, new_y)

Now you have two-dimensional points. Run the code above to sort them.

Answer 2

After many hours, I finally found a solution. @Pranav Hosangadi's solution worked for the 2D side of things. However, I was having trouble projecting the 3D coordinates to 2D coordinates using the second part of his solution. I also tried projecting the coordinates as described in this answer, but it did not work as intended. I then discovered an API function called location_3d_to_region_2d() (see docs) which, as the name implies, gets the 2D screen coordinates in pixels of the given 3D coordinate. I didn't need to necessarily "project" anything into 2D in the first place, getting the screen coordinates worked perfectly fine and is much more simple. From that point, I could sort the coordinates using Pranav's function with some slight adjustments to get it in the order illustrated in the screenshot of my first post and I wanted it returned as a list instead of a NumPy array.

import bpy
from bpy_extras.view3d_utils import location_3d_to_region_2d
import numpy

def sort_points(pts):
    """Sort 4 points in a winding order"""
    pts = numpy.array(pts)
    centroid = numpy.sum(pts, axis=0) / pts.shape[0]
    vector_from_centroid = pts - centroid
    vector_angle = numpy.arctan2(
        vector_from_centroid[:, 1], vector_from_centroid[:, 0])
    # Find the indices that give a sorted vector_angle array
    sort_order = numpy.argsort(-vector_angle)

    # Apply sort_order to original pts array.
    return list(sort_order)

# Get 2D screen coords of selected vertices
region = bpy.context.region
region_3d = bpy.context.space_data.region_3d

corners2d = []
for corner in selected_verts:
    corners2d.append(location_3d_to_region_2d(
        region, region_3d, corner))

# Sort the 2d points in a winding order
sort_order = sort_points(corners2d)
sorted_corners = [selected_verts[i] for i in sort_order]

Thanks, Pranav for your time and patience in helping me solve this problem!

Answer 3

There is a simpler and faster solution for the Blender case:

1.) The following code sorts 4 planar points in 2D (vertices of the plane object in Blender) very efficiently:

def sort_clockwise(pts):
    rect = np.zeros((4, 2), dtype="float32")
    s = pts.sum(axis=1)
    rect[0] = pts[np.argmin(s)]
    rect[2] = pts[np.argmax(s)]
    diff = np.diff(pts, axis=1)
    rect[1] = pts[np.argmin(diff)]
    rect[3] = pts[np.argmax(diff)]
    return rect

2.) Blender keeps vertices related data, such as the translation, rotation and scale in the world matrix. If you query for vertices.co(ordinates) only, you just get the original coordinates, without translation, rotation and scaling. But that does not affect the order of vertices. That simplifies the problem because what you get is actually a 2D (with z's = 0) mesh data. If you sort that 2D data (excluding z's) you will get the information, the sort indices for the 3D sorted data. You can modify the code above to get the indices from that 2D array. For the plane object of Blender, for some reason the order is always [0,1,3,2], not [0,1,2,3]. The following modified code gives the sorted indices for the vertices data in 2D.

def sorted_ix_clockwise(pts):
    #rect = zeros((4, 2), dtype="float32")
    ix = array([0,0,0,0])
    s = pts.sum(axis=1)
    #rect[0] = pts[argmin(s)]
    #rect[2] = pts[argmax(s)]
    ix[0] = argmin(s)
    ix[2] = argmax(s)
    dif = diff(pts, axis=1)
    #rect[1] = pts[argmin(dif)]
    #rect[3] = pts[argmax(dif)]
    ix[1] = argmin(dif)
    ix[3] = argmax(dif)
    return ix

You can use these indices to get the actual 3D sorted data, which you can obtain by multiplying vertices coordinates with the world matrix to include any translation, rotation and scaling.