Sorting 3D points on the surface of a sphere in clockwise order

Question

I have an array of 3D points on the surface of the unit sphere, and a center point C (also on the surface of the unit sphere). How can I sort these points so they are in clockwise order relative to the surface?

Answer 1

This question may have been better placed on http://math.stackexchange.com or https://mathoverflow.net/ since it relates more to linear algebra than how to code it.

At any rate, here is what I did for ordering points into convex spherical polygons. There will probably be cases with non-convex polygons where it will break down, but I haven't checked as I am working with convex polygons exclusively for what I needed to do. So, an incomplete answer but hopefully a useful one for somebody else. This give a better solution than what Good Luck suggested since it can deal with points around a pole or crossing the antimeridian.

A few assumptions: I am working on a unit sphere, where all points of interest can be thought of as vectors from the origin and also represent vertices of a convex spherical polygon. Also, since the set of points can describe two potential spherical polygons, I am assuming that the smaller one is the one we are after. There are N points/vertices of interest held in an Nx3 array (so each vector v[0,:], v[1,:], ... v[N-1,:] is defined with (x,y,z) coordinates and has a norm of 1. In relation to user3235832's question, I am assuming the axis of interest is somewhere inside the polygon.

Step 1: Pick a location v_c that is somewhere inside the polygon. I simply normalized the mean value of my points. The normalization ensures it too is on the sphere. For non-convex polygons there may be a need to do something more robust to ensure the point would be inside the polygon.

Step 2: Arbitrarily using the first point as an anchor, find the interior angle between that point, our centre position, and each of the other vertices. The interior angle is also the angle between the two planes containing the great circles for each edge. So that would be something like

import numpy as np
alpha = np.zeros(N-1) # initialize array
for i in range(1,N):
   alpha[i-1] = np.arccos(np.dot(np.cross(v[0], v_c),np.cross(v_c, v[i])))

However, this may give the wrong value for the arccos (e.g. pi/3 instead of 2pi - pi/3) so it is necessary to check whether we've got the correct angle here. To do this, I check if the angle is clockwise or counterclockwise using the approach suggested here:

import numpy as np
def is_ccw(v_0,v_c,v_i):
    # checks if the smaller interior angle for the great circles connecting u-v and v-w is CCW
    return(np.dot(np.cross(v_c - v_0,v_i - v_c), v_i) < 0)

If the the angle is not CCW then it would need to be corrected, i.e. alpha = 2 * np.pi - alpha.

Step 3: Sort the points by angle from v_1.

Answer 2

A bit late but you can calculate the normal and then see if along the normal axis they average clockwise or counter clockwise around the normal or around any single point of the surface. The second way you have one less point to track but it's basically the same.

1. calculate the center point of the surface (the sum by coordinate divided by the count by coordinates);
2. subtract the center point from the surface points;
3. rotate the surface points by the surface normal; ** you are only dealing with 2 axis now because z will always be zero;
4. sum the relative angles (sum(atan(p[i] - p[0]) for i in (1,2)))
5. if it's positive then it's clockwise.

some code...

from math import sqrt, atan2, pi
twopi = 2 * pi
dists = lambda s: sqrt(sum(c*c for c in s))

import numpy as np
trans_y_mat = lambda dx, dz: np.array(((dz,0,dx),(0,1,0),(-dx,0,dz)) \
    , dtype=np.float64)


def p3d_normal(p3d):
  Ux, Uy, Uz = (p3d[1][i] - p3d[0][i] for i in (0,1,2))
  Vx, Vy, Vz = (p3d[2][i] - p3d[0][i] for i in (0,1,2))
  N = (Uy*Vz - Uz*Vy, Uz*Vx - Ux*Vz, Ux*Vy - Uy*Vx)
  d = dists(N)
  return tuple(c / d for c in N)


def p3d_is_clockwise(p3d, p3n=None):
  if p3n is None: p3n = p3d_normal(p3d)
  dnx, dnz = p3n[0]/p3n[1], p3n[2]/p3n[1]
  dn = dists((dnz, dnx))
  mn = trans_y_mat(dnz/dn, dnx/dn)
  p2d = np.matmul(mn, p3d)[:,:2]
  asum = 0.0
  xp,yp = p2d[0]
  ap = 0.0
  for (xn,yn) in p2d[1:]:
    an = atan2(yn-yp, xn-xp)
    asum += (an - ap + pi) % twopi - pi
    xp, yp, an = xn, yp, ap
  return asum >= 0


faces = (((0.0, 0.0, 100.0), (30.9017, 0.0, 95.1057), (15.4508, 26.7617, 95.1057)) \
    , ((0.0, 0.0, 100.0), (15.4508, 26.7617, 95.1057), (-15.4508, 26.7617, 95.1057)) \
    , ((0.0, 0.0, 100.0), (-15.4508, 26.7617, 95.1057), (-30.9017, 0.0, 95.1057)) \
    , ((0.0, 0.0, 100.0), (-30.9017, 0.0, 95.1057), (-15.4508, -26.7617, 95.1057)) \
    , ((0.0, 0.0, 100.0), (-15.4508, -26.7617, 95.1057), (15.4508, -26.7617, 95.1057)) \
    , ((0.0, 0.0, 100.0), (15.4508, -26.7617, 95.1057), (30.9017, -0.0, 95.1057)) \
    , ((30.9017, 0.0, 95.1057), (50.9037, 29.3893, 80.9017), (58.7785, 0.0, 80.9017)) \
    , ((30.9017, 0.0, 95.1057), (15.4508, 26.7617, 95.1057), (50.9037, 29.3893, 80.9017)) \
    , ((15.4508, 26.7617, 95.1057), (29.3893, 50.9037, 80.9017), (50.9037, 29.3893, 80.9017)) \
    , ((15.4508, 26.7617, 95.1057), (0.0, 58.7785, 80.9017), (29.3893, 50.9037, 80.9017)))

for face in faces:
  print(p3d_is_clockwise(face))

Output:

False
False
True
False
False
False
False
False
True
True

Answer 3

Another very late answer, I found this post as I had the same problem and found myself using this form here to rubber duck it out.

I have a solve that is very specific to my needs, so it takes some shortcuts for the sake of speed that make it less general purpose. I'm recursively subdividing an icosahedron into a geosphere, each point has a list of neighbours I need to reorder clockwise so that the polygons face the right way when rendered. I figured this task is the most likely reason anyone would need to clockwise order points on a unit sphere, so I might as well post the method I came up with.

So in my specific case, each point has 6 neighbours that are fairly evenly spaced, except for 12 points from the original icosahedron that have 5 evenly spaced neighbours. It doesn't matter which of the neighbours is first in the sorted list, because the last will attach to the first anyway, only the clockwise order matters.

I have one centre point C, 5 or 6 neighbours N[ ], and an output array S[ ]. I subtract C from all neighbours to get deltas around an origin and normalise those into D[ ]. D[0] becomes a Y axis of sorts, and the cross product of C and D[0] becomes the X axis.

Because there is always one and only one point in each quadrant, I can take the dot product of each point D[1 thru 5] vs Y and X and based on whether those are positive or negative, determine how many indices to offset each point from N[0] and in which direction to offset. So I end up putting N[0] straight into S[2] then looping through N setting S[2+offset*direction] = N. where offset is 1 or 2 and direction is -1 or 1.

This handles 5 of the points, but when there are 6 points, one will be close to directly opposite N[0] so the Dot product vs X is near 0, and vs Y it's near -1, I test for that manually and force the offset to +3 in that case.

Hopefully this might help the next person who goes down this path, it's maybe not the best way, maybe a bit of a kludge, but reasonably efficient as it tests one vs. many instead of testing many vs. many by making a lot of assumptions about the input.