def project_3d_point(p: Vector):
    region_data = bpy.context.region_data

    # Get the two components to calculate M
    modelview_matrix = region_data.view_matrix
    projection_matrix = region_data.perspective_matrix

    # Compute P’ = M * P
    p1 = projection_matrix @ modelview_matrix @ Vector((p.x, p.y, p.z, 1))

    # Normalize in: x’’ = x’ / w’, y’’ = y’ / w’
    p2 = Vector(((p1.x/p1.w, p1.y/p1.w)))

    return p2

# Project points onto 2D plane
corners2d = []
for vert in selected_verts:
    corners2d.append(project_3d_point(vert))

# Sort the 2d points in a winding order
sort_order = list(sort_points(np.array(corners2d)))
sorted_corners = [None] * 4
sorted_corners[0] = selected_verts[sort_order.index(0)]
sorted_corners[1] = selected_verts[sort_order.index(1)]
sorted_corners[2] = selected_verts[sort_order.index(2)]
sorted_corners[3] = selected_verts[sort_order.index(3)]