定位球形多边形的质心（质心）

我正在试图找出如何最好地定位覆盖在单位球面上的任意形状的质心，input顺序（顺时针或反cw）顶点的形状边界。 顶点的密度沿着边界是不规则的，所以它们之间的弧长一般是不相等的。 由于形状可能非常大（半个半球），所以通常不可能简单地将顶点投影到平面上，并使用平面方法，详见维基百科（对不起，我不允许超过2个超链接作为新手）。 稍微好一点的方法是使用在球坐标系中操作的平面几何graphics，但是如果使用大的多边形，这种方法也会失败，正如这里很好地说明的那样。 在同一页上，“Cffk”突出介绍了一种计算球形三angular形质心的方法。 我试图实现这个方法，但没有成功，我希望有人可以发现这个问题？

我已经将variables的定义与文件中的类似，以便于比较。 input（数据）是经度/纬度坐标列表，由代码转换为[x，y，z]坐标。 对于每个三angular形，我已经任意固定一个点作为+ z极点，其他两个顶点由沿着多边形边界的一对相邻点组成。 代码沿边界步进（从任意点开始），依次使用多边形的每个边界段作为三angular形边。 为每个这些单独的球形三angular形确定一个子质心，并根据三angular形面积对它们进行加权，并将其相加以计算总的多边形质心。 运行代码时我没有遇到任何错误，但返回的总质心显然是错误的（我已经运行了质心位置明确的一些非常基本的形状）。 我还没有find任何明智的模式，在重心的位置返回…所以目前我不知道什么是错误的，无论是在math或代码（虽然，怀疑是math）。

下面的代码应该工作原样复制粘贴，如果你想尝试。 如果你已经安装了matplotlib和numpy，它会绘制结果（如果你不这样做，它将忽略绘图）。 您只需要将代码下方的经度/纬度数据放入一个名为example.txt的文本文件中即可。

from math import * try: import matplotlib as mpl import matplotlib.pyplot from mpl_toolkits.mplot3d import Axes3D import numpy plotting_enabled = True except ImportError: plotting_enabled = False def sph_car(point): if len(point) == 2: point.append(1.0) rlon = radians(float(point[0])) rlat = radians(float(point[1])) x = cos(rlat) * cos(rlon) * point[2] y = cos(rlat) * sin(rlon) * point[2] z = sin(rlat) * point[2] return [x, y, z] def xprod(v1, v2): x = v1[1] * v2[2] - v1[2] * v2[1] y = v1[2] * v2[0] - v1[0] * v2[2] z = v1[0] * v2[1] - v1[1] * v2[0] return [x, y, z] def dprod(v1, v2): dot = 0 for i in range(3): dot += v1[i] * v2[i] return dot def plot(poly_xyz, g_xyz): fig = mpl.pyplot.figure() ax = fig.add_subplot(111, projection='3d') # plot the unit sphere u = numpy.linspace(0, 2 * numpy.pi, 100) v = numpy.linspace(-1 * numpy.pi / 2, numpy.pi / 2, 100) x = numpy.outer(numpy.cos(u), numpy.sin(v)) y = numpy.outer(numpy.sin(u), numpy.sin(v)) z = numpy.outer(numpy.ones(numpy.size(u)), numpy.cos(v)) ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', linewidth=0, alpha=0.3) # plot 3d and flattened polygon x, y, z = zip(*poly_xyz) ax.plot(x, y, z) ax.plot(x, y, zs=0) # plot the alleged 3d and flattened centroid x, y, z = g_xyz ax.scatter(x, y, z, c='r') ax.scatter(x, y, 0, c='r') # display ax.set_xlim3d(-1, 1) ax.set_ylim3d(-1, 1) ax.set_zlim3d(0, 1) mpl.pyplot.show() lons, lats, v = list(), list(), list() # put the two-column data at the bottom of the question into a file called # example.txt in the same directory as this script with open('example.txt') as f: for line in f.readlines(): sep = line.split() lons.append(float(sep[0])) lats.append(float(sep[1])) # convert spherical coordinates to cartesian for lon, lat in zip(lons, lats): v.append(sph_car([lon, lat, 1.0])) # z unit vector/pole ('north pole'). This is an arbitrary point selected to act as one #(fixed) vertex of the summed spherical triangles. The other two vertices of any #triangle are composed of neighboring vertices from the polygon boundary. np = [0.0, 0.0, 1.0] # Gx,Gy,Gz are the cartesian coordinates of the calculated centroid Gx, Gy, Gz = 0.0, 0.0, 0.0 for i in range(-1, len(v) - 1): # cycle through the boundary vertices of the polygon, from 0 to n if all((v[i][0] != v[i+1][0], v[i][1] != v[i+1][1], v[i][2] != v[i+1][2])): # this just ignores redundant points which are common in my larger input files # A,B,C are the internal angles in the triangle: 'np-v[i]-v[i+1]-np' A = asin(sqrt((dprod(np, xprod(v[i], v[i+1])))**2 / ((1 - (dprod(v[i+1], np))**2) * (1 - (dprod(np, v[i]))**2)))) B = asin(sqrt((dprod(v[i], xprod(v[i+1], np)))**2 / ((1 - (dprod(np , v[i]))**2) * (1 - (dprod(v[i], v[i+1]))**2)))) C = asin(sqrt((dprod(v[i + 1], xprod(np, v[i])))**2 / ((1 - (dprod(v[i], v[i+1]))**2) * (1 - (dprod(v[i+1], np))**2)))) # A/B/Cbar are the vertex angles, such that if 'O' is the sphere center, Abar # is the angle (v[i]-Ov[i+1]) Abar = acos(dprod(v[i], v[i+1])) Bbar = acos(dprod(v[i+1], np)) Cbar = acos(dprod(np, v[i])) # e is the 'spherical excess', as defined on wikipedia e = A + B + C - pi # mag1/2/3 are the magnitudes of vectors np,v[i] and v[i+1]. mag1 = 1.0 mag2 = float(sqrt(v[i][0]**2 + v[i][1]**2 + v[i][2]**2)) mag3 = float(sqrt(v[i+1][0]**2 + v[i+1][1]**2 + v[i+1][2]**2)) # vec1/2/3 are cross products, defined here to simplify the equation below. vec1 = xprod(np, v[i]) vec2 = xprod(v[i], v[i+1]) vec3 = xprod(v[i+1], np) # multiplying vec1/2/3 by e and respective internal angles, according to the #posted paper for x in range(3): vec1[x] *= Cbar / (2 * e * mag1 * mag2 * sqrt(1 - (dprod(np, v[i])**2))) vec2[x] *= Abar / (2 * e * mag2 * mag3 * sqrt(1 - (dprod(v[i], v[i+1])**2))) vec3[x] *= Bbar / (2 * e * mag3 * mag1 * sqrt(1 - (dprod(v[i+1], np)**2))) Gx += vec1[0] + vec2[0] + vec3[0] Gy += vec1[1] + vec2[1] + vec3[1] Gz += vec1[2] + vec2[2] + vec3[2] approx_expected_Gxyz = (0.78, -0.56, 0.27) print('Approximate Expected Gxyz: {0}\n' ' Actual Gxyz: {1}' ''.format(approx_expected_Gxyz, (Gx, Gy, Gz))) if plotting_enabled: plot(v, (Gx, Gy, Gz)) 

在此先感谢您的任何build议或见解。

编辑：这里是一个图，显示了一个多边形的单位球面的投影和我从代码计算的结果质心。 显然，质心是错误的，因为多边形相当小而凸，但质心落在其外围。

编辑：这是一个非常类似的一套坐标，但在原来的[lon，lat]格式我通常使用（现在被更新的代码转换为[x，y，z]）。

  -39.366295 -1.633460 -47.282630 -0.740433 -53.912136 0.741380 -59.004217 2.759183 -63.489005 5.426812 -68.566001 8.712068 -71.394853 11.659135 -66.629580 15.362600 -67.632276 16.827507 -66.459524 19.069327 -63.819523 21.446736 -61.672712 23.532143 -57.538431 25.947815 -52.519889 28.691766 -48.606227 30.646295 -45.000447 31.089437 -41.549866 32.139873 -36.605156 32.956277 -32.010080 34.156692 -29.730629 33.756566 -26.158767 33.714080 -25.821513 34.179648 -23.614658 36.173719 -20.896869 36.977645 -17.991994 35.600074 -13.375742 32.581447 -9.554027 28.675497 -7.825604 26.535234 -7.825604 26.535234 -9.094304 23.363132 -9.564002 22.527385 -9.713885 22.217165 -9.948596 20.367878 -10.496531 16.486580 -11.151919 12.666850 -12.350144 8.800367 -15.446347 4.993373 -20.366139 1.132118 -24.784805 -0.927448 -31.532135 -1.910227 -39.366295 -1.633460 

编辑：更多的例子…有4个顶点，定义一个完美的平方[1,0,0]我得到预期的结果： 然而，从一个不对称的三angular形中，我得到的是一个不重要的质心…质心实际上落在球体的远侧（这里作为对映体投射到前侧）： 有趣的是，如果我反转列表（从顺时针到逆时针顺序，反之亦然），则质心估计看起来是“稳定的”，质心相应地反转。

Answer 1

我认为这将做到这一点。 你应该能够通过复制下面的代码复制这个结果。

• 您需要将纬度和经度数据放在一个名为longitude and latitude.txt的文件中。 您可以复制粘贴代码下面包含的原始示例数据。
• 如果你有mplotlib它会另外产生下面的情节
• 对于不明显的计算，我包括一个链接，解释发生了什么
• 在下面的图表中，参考vector非常短（r = 1/10），因此3d-质心更容易看清。 您可以轻松删除缩放以最大化准确性。
• 注意op：我重写了几乎所有的东西，所以我不确定原始代码不能正确工作。 不过，至less我认为没有考虑到处理顺时针/逆时针三angular形顶点的需要。

---

传说：

• （黑线）参考vector
• （小红点）球形三angular形3d-质心
• （大红色/蓝色/绿色点）3d质心/投影到表面/投影到xy平面
• （蓝/绿线）球面多边形和投影到xy平面上

---

 from math import * try: import matplotlib as mpl import matplotlib.pyplot from mpl_toolkits.mplot3d import Axes3D import numpy plotting_enabled = True except ImportError: plotting_enabled = False def main(): # get base polygon data based on unit sphere r = 1.0 polygon = get_cartesian_polygon_data(r) point_count = len(polygon) reference = ok_reference_for_polygon(polygon) # decompose the polygon into triangles and record each area and 3d centroid areas, subcentroids = list(), list() for ia, a in enumerate(polygon): # build an abc point set ib = (ia + 1) % point_count b, c = polygon[ib], reference if points_are_equivalent(a, b, 0.001): continue # skip nearly identical points # store the area and 3d centroid areas.append(area_of_spherical_triangle(r, a, b, c)) tx, ty, tz = zip(a, b, c) subcentroids.append((sum(tx)/3.0, sum(ty)/3.0, sum(tz)/3.0)) # combine all the centroids, weighted by their areas total_area = sum(areas) subxs, subys, subzs = zip(*subcentroids) _3d_centroid = (sum(a*subx for a, subx in zip(areas, subxs))/total_area, sum(a*suby for a, suby in zip(areas, subys))/total_area, sum(a*subz for a, subz in zip(areas, subzs))/total_area) # shift the final centroid to the surface surface_centroid = scale_v(1.0 / mag(_3d_centroid), _3d_centroid) plot(polygon, reference, _3d_centroid, surface_centroid, subcentroids) def get_cartesian_polygon_data(fixed_radius): cartesians = list() with open('longitude and latitude.txt') as f: for line in f.readlines(): spherical_point = [float(v) for v in line.split()] if len(spherical_point) == 2: spherical_point.append(fixed_radius) cartesians.append(degree_spherical_to_cartesian(spherical_point)) return cartesians def ok_reference_for_polygon(polygon): point_count = len(polygon) # fix the average of all vectors to minimize float skew polyx, polyy, polyz = zip(*polygon) # /10 is for visualization. Remove it to maximize accuracy return (sum(polyx)/(point_count*10.0), sum(polyy)/(point_count*10.0), sum(polyz)/(point_count*10.0)) def points_are_equivalent(a, b, vague_tolerance): # vague tolerance is something like a percentage tolerance (1% = 0.01) (ax, ay, az), (bx, by, bz) = a, b return all(((ax-bx)/ax < vague_tolerance, (ay-by)/ay < vague_tolerance, (az-bz)/az < vague_tolerance)) def degree_spherical_to_cartesian(point): rad_lon, rad_lat, r = radians(point[0]), radians(point[1]), point[2] x = r * cos(rad_lat) * cos(rad_lon) y = r * cos(rad_lat) * sin(rad_lon) z = r * sin(rad_lat) return x, y, z def area_of_spherical_triangle(r, a, b, c): # points abc # build an angle set: A(CAB), B(ABC), C(BCA) # http://math.stackexchange.com/a/66731/25581 A, B, C = surface_points_to_surface_radians(a, b, c) E = A + B + C - pi # E is called the spherical excess area = r**2 * E # add or subtract area based on clockwise-ness of abc # http://stackoverflow.com/a/10032657/377366 if clockwise_or_counter(a, b, c) == 'counter': area *= -1.0 return area def surface_points_to_surface_radians(a, b, c): """build an angle set: A(cab), B(abc), C(bca)""" points = a, b, c angles = list() for i, mid in enumerate(points): start, end = points[(i - 1) % 3], points[(i + 1) % 3] x_startmid, x_endmid = xprod(start, mid), xprod(end, mid) ratio = (dprod(x_startmid, x_endmid) / ((mag(x_startmid) * mag(x_endmid)))) angles.append(acos(ratio)) return angles def clockwise_or_counter(a, b, c): ab = diff_cartesians(b, a) bc = diff_cartesians(c, b) x = xprod(ab, bc) if x < 0: return 'clockwise' elif x > 0: return 'counter' else: raise RuntimeError('The reference point is in the polygon.') def diff_cartesians(positive, negative): return tuple(p - n for p, n in zip(positive, negative)) def xprod(v1, v2): x = v1[1] * v2[2] - v1[2] * v2[1] y = v1[2] * v2[0] - v1[0] * v2[2] z = v1[0] * v2[1] - v1[1] * v2[0] return [x, y, z] def dprod(v1, v2): dot = 0 for i in range(3): dot += v1[i] * v2[i] return dot def mag(v1): return sqrt(v1[0]**2 + v1[1]**2 + v1[2]**2) def scale_v(scalar, v): return tuple(scalar * vi for vi in v) def plot(polygon, reference, _3d_centroid, surface_centroid, subcentroids): fig = mpl.pyplot.figure() ax = fig.add_subplot(111, projection='3d') # plot the unit sphere u = numpy.linspace(0, 2 * numpy.pi, 100) v = numpy.linspace(-1 * numpy.pi / 2, numpy.pi / 2, 100) x = numpy.outer(numpy.cos(u), numpy.sin(v)) y = numpy.outer(numpy.sin(u), numpy.sin(v)) z = numpy.outer(numpy.ones(numpy.size(u)), numpy.cos(v)) ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', linewidth=0, alpha=0.3) # plot 3d and flattened polygon x, y, z = zip(*polygon) ax.plot(x, y, z, c='b') ax.plot(x, y, zs=0, c='g') # plot the 3d centroid x, y, z = _3d_centroid ax.scatter(x, y, z, c='r', s=20) # plot the spherical surface centroid and flattened centroid x, y, z = surface_centroid ax.scatter(x, y, z, c='b', s=20) ax.scatter(x, y, 0, c='g', s=20) # plot the full set of triangular centroids x, y, z = zip(*subcentroids) ax.scatter(x, y, z, c='r', s=4) # plot the reference vector used to findsub centroids x, y, z = reference ax.plot((0, x), (0, y), (0, z), c='k') ax.scatter(x, y, z, c='k', marker='^') # display ax.set_xlim3d(-1, 1) ax.set_ylim3d(-1, 1) ax.set_zlim3d(0, 1) mpl.pyplot.show() # run it in a function so the main code can appear at the top main() 

---

这里是可以粘贴到longitude and latitude.txt经纬度数据

  -39.366295 -1.633460 -47.282630 -0.740433 -53.912136 0.741380 -59.004217 2.759183 -63.489005 5.426812 -68.566001 8.712068 -71.394853 11.659135 -66.629580 15.362600 -67.632276 16.827507 -66.459524 19.069327 -63.819523 21.446736 -61.672712 23.532143 -57.538431 25.947815 -52.519889 28.691766 -48.606227 30.646295 -45.000447 31.089437 -41.549866 32.139873 -36.605156 32.956277 -32.010080 34.156692 -29.730629 33.756566 -26.158767 33.714080 -25.821513 34.179648 -23.614658 36.173719 -20.896869 36.977645 -17.991994 35.600074 -13.375742 32.581447 -9.554027 28.675497 -7.825604 26.535234 -7.825604 26.535234 -9.094304 23.363132 -9.564002 22.527385 -9.713885 22.217165 -9.948596 20.367878 -10.496531 16.486580 -11.151919 12.666850 -12.350144 8.800367 -15.446347 4.993373 -20.366139 1.132118 -24.784805 -0.927448 -31.532135 -1.910227 -39.366295 -1.633460 

Answer 2

我认为一个好的近似值是使用加权笛卡尔坐标计算质心，并将结果投影到球面上（假设坐标原点是(0, 0, 0)^T ）。

令(p[0], p[1], ... p[n-1])多边形的n个点。 近似（笛卡尔）质心可以通过以下公式计算：

 c = 1 / w * (sum of w[i] * p[i]) 

而w是所有权重的总和，而p[i]是多边形点， w[i]是该点的权重，例如

 w[i] = |p[i] - p[(i - 1 + n) % n]| / 2 + |p[i] - p[(i + 1) % n]| / 2 

而|x| 是向量x的长度。 也就是说，一个点的权重是前一半长度和下一个多边形点长度的一半。

这个质心c现在可以通过以下方式投影到球体上：

 c' = r * c / |c| 

而r是球体的半径。

要考虑多边形的方向(ccw, cw) ，结果可能是

 c' = - r * c / |c|. 

Answer 3

为了澄清：感兴趣的数量是真实的三维质心（即三维质心，即三维中心）在单位球面上的投影。

由于你所关心的是从原点到三维质心的方向，所以你根本不需要打扰区域; 计算时刻（即三维重心时间面积）更容易。 单位球体上闭合path左侧的区域的矩是当你在path中行走时左边单位vector的一半。 这是由斯托克斯定理的一个非显而易见的应用引起的。 见http://www.owlnet.rice.edu/~fjones/chap13.pdf问题13-12。;

特别是，对于一个球形多边形，矩是（axb）/ || axb ||的和的一半 *（a和b之间的angular度）为每对连续顶点a，b。 （这是path左侧的区域;否则为path右侧的区域）。

（如果你真的想要三维质心，只需要计算面积并将其分成两部分，比较面积也可能有助于select哪个区域称为“多边形”。）

这是一些代码; 这很简单：

 #!/usr/bin/python import math def plus(a,b): return [x+y for x,y in zip(a,b)] def minus(a,b): return [xy for x,y in zip(a,b)] def cross(a,b): return [a[1]*b[2]-a[2]*b[1], a[2]*b[0]-a[0]*b[2], a[0]*b[1]-a[1]*b[0]] def dot(a,b): return sum([x*y for x,y in zip(a,b)]) def length(v): return math.sqrt(dot(v,v)) def normalized(v): l = length(v); return [1,0,0] if l==0 else [x/l for x in v] def addVectorTimesScalar(accumulator, vector, scalar): for i in xrange(len(accumulator)): accumulator[i] += vector[i] * scalar def angleBetweenUnitVectors(a,b): # http://www.plunk.org/~hatch/rightway.php if dot(a,b) < 0: return math.pi - 2*math.asin(length(plus(a,b))/2.) else: return 2*math.asin(length(minus(a,b))/2.) def sphericalPolygonMoment(verts): moment = [0.,0.,0.] for i in xrange(len(verts)): a = verts[i] b = verts[(i+1)%len(verts)] addVectorTimesScalar(moment, normalized(cross(a,b)), angleBetweenUnitVectors(a,b) / 2.) return moment if __name__ == '__main__': import sys def lonlat_degrees_to_xyz(lon_degrees,lat_degrees): lon = lon_degrees*(math.pi/180) lat = lat_degrees*(math.pi/180) coslat = math.cos(lat) return [coslat*math.cos(lon), coslat*math.sin(lon), math.sin(lat)] verts = [lonlat_degrees_to_xyz(*[float(v) for v in line.split()]) for line in sys.stdin.readlines()] #print "verts = "+`verts` moment = sphericalPolygonMoment(verts) print "moment = "+`moment` print "centroid unit direction = "+`normalized(moment)` 

对于示例多边形，这给出了答案（单位vector）：

 [-0.7644875430808217, 0.579935445918147, -0.2814847687566214] 

这与@ KobeJohn的代码所计算的答案大致相同，但是更准确，该代码使用粗略的容差和平面逼近子中心：

 [0.7628095787179151, -0.5977153368303585, 0.24669398601094406] 

这两个答案的方向大致相反（所以我猜在这种情况下，科比约翰的代码决定把这个区域放在path的右边 ）。