按顺时针顺序对球体表面上的 3D 点进行排序

Question

我在单位球体的表面上有一个 3D 点数组，还有一个中心点 C（也在单位球体的表面上）。如何对这些点进行排序，使它们相对于表面按顺时针顺序排列？

Answer 1

这个问题最好放在http://math.stackexchange.com 或https://mathoverflow.net/ 上，因为它更多地与线性代数相关，而不是如何编码。

无论如何，这就是我为将点排序为凸球面多边形所做的工作。可能存在非凸多边形会分解的情况，但我没有检查过，因为我正在使用凸多边形专门用于我需要做的事情。所以，一个不完整的答案，但希望对其他人有用。这提供了比 Good Luck 建议的更好的解决方案，因为它可以处理极点周围的点或穿过反子午线。

一些假设：我正在研究一个单位球体，其中所有兴趣点都可以被认为是来自原点的向量，也可以表示凸球面多边形的顶点。此外，由于这组点可以描述两个潜在的球形多边形，我假设较小的那个是我们所追求的那个。 Nx3 数组中有 N 个感兴趣的点/顶点（因此每个向量 v[0,:], v[1,:], ... v[N-1,:] 用 (x,y,z) 坐标定义，范数为 1。关于 user3235832 的问题，我假设感兴趣的轴在多边形内的某个地方。

第 1 步：选择多边形内某处的位置 v_c。我只是将我的分数的平均值标准化。标准化确保它也在球体上。对于非凸多边形，可能需要做一些更稳健的事情来确保该点位于多边形内部。

第 2 步： 任意使用第一个点作为锚点，找到该点、我们的中心位置和其他每个顶点之间的内角。内角也是包含每个边的大圆的两个平面之间的角度。所以会是这样的

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

但是，这可能会为 arccos 提供错误的值（例如 pi/3 而不是 2pi - pi/3），因此有必要检查我们是否在这里获得了正确的角度。为此，我使用建议的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)

如果角度不是逆时针，则需要进行校正，即alpha = 2 * np.pi - alpha。

第 3 步：从 v_1 按角度对点进行排序。

Answer 2

有点晚了，但是您可以计算法线，然后查看沿法线轴的平均值是顺时针还是逆时针围绕法线或曲面的任何单个点。第二种方法你要跟踪的点少了，但基本上是一样的。

1. 计算曲面的中心点（坐标总和除以坐标计数）；
2. 从表面点中减去中心点；
3. 按表面法线旋转表面点； ** 您现在只处理 2 轴，因为 z 将始终为零；
4. 对相对角度求和 (sum(atan(p[i] - p[0]) for i in (1,2)))
5. 如果为正，则为顺时针。

一些代码...

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

输出：

False
False
True
False
False
False
False
False
True
True