sort points

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])
    # Find the indices that give a sorted vector_angle array
    sort_order = np.argsort(vector_angle)

    # Apply sort_order to original pts array.
    return pts[sort_order, :]

# selected_verts are the points that the user selects...

view_mat = bpy.context.region_data.view_matrix
# Project points onto 2D plane
corners = np.array(
    [view_mat @ p for p in selected_verts])
plane_normal = np.cross(
    corners[1, :] - corners[0, :], corners[2, :] - corners[0, :])
plane_normal = plane_normal / np.linalg.norm(plane_normal)

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

new_y = np.cross(plane_normal, new_x)

# Maps the indices of the original 3d points to the projected 2d points
corners2dmap = {}
for index, corner in enumerate(corners):
    proj_x = Vector.dot(
        Vector(corner) - Vector(new_origin), Vector(new_x))
    proj_y = Vector.dot(
        Vector(corner) - Vector(new_origin), Vector(new_y))
    # Cannot use Blender's Vector type as a key for a dict so must convert to tuple for the mapping
    corners2dmap[tuple((proj_x, proj_y))] = index

# Sort the 2d points in a winding order
sorted_corners2d = sort_points(np.array(list(corners2dmap.keys())))
sorted_corners = []
# Obtain the sorted list of 3d points with the view_matrix transformation reversed
for corner in sorted_corners2d:
    sorted_corners.append(view_mat.inverted() @
                            Vector(corners[corners2dmap[tuple(corner)]]))