Localisation du centre de masse des polygones sphériques

Question

Je suis en train de travailler sur la meilleure façon de localiser le centre de gravité d'une forme arbitraire drapé sur une sphère unité, avec l'entrée en cours de commande (dans le sens horaire ou anti-cw) de sommets pour les limites de la forme. La densité de sommets est irrégulière le long de la frontière, de sorte que l'arc-longueurs entre eux ne sont généralement pas égale. Parce que les formes peuvent être très grandes (la moitié d'un hémisphère) il n'est généralement pas possible de simplement de projet les sommets d'un avion de ligne et l'utilisation plane méthodes, comme détaillé sur Wikipédia (désolé, je n'ai pas le droit à plus de 2 liens hypertexte en tant que nouveau venu). Une légèrement meilleure approche implique l'utilisation de la géométrie plane manipulés en coordonnées sphériques, mais encore une fois, avec de grands polygones cette méthode échoue, joliment illustré ici. Sur la même page, 'Cffk a mis en lumière le présent document qui décrit une méthode pour le calcul du centre de gravité des triangles sphériques. J'ai essayé de mettre en œuvre cette méthode, mais sans succès, et j'espère que quelqu'un peut repérer le problème?

J'ai gardé les définitions de variables similaires à ceux dans le papier pour le rendre plus facile à comparer. L'entrée (données) est une liste de coordonnées latitude/longitude coordonnées, converti à [x,y,z] coordonnées par le code. Pour chacun des triangles, j'ai fixé arbitrairement un point pour être le +z-pole, les deux autres sommets étant composé d'une paire de voisins points le long de la limite de polygone. Le code suit le long de la frontière (en commençant à un point arbitraire), à l'aide de chaque segment limite du polygone comme un triangle de côté à son tour. Un sous-centre de gravité est déterminé pour chacune de ces triangles sphériques et ils sont pondérés en fonction de l'aire du triangle et ajouté pour calculer le total polygone centre de gravité. Je n'ai pas d'erreur lors de l'exécution du code, mais le total des centroïdes retournés sont clairement mal (j'ai des formes là où le centre de gravité de l'emplacement est sans ambiguïté). Je n'ai trouvé aucun motif raisonnable dans l'emplacement des centroïdes retourné...donc pour le moment je ne suis pas sûr de ce qui ne va pas, que ce soit en maths ou en code (bien que, la suspicion est les maths).

Le code ci-dessous devrait fonctionner copier-coller comme c'est si vous voulez l'essayer. Si vous avez matplotlib et numpy installé, il va tracer les résultats (il ignore le traçage si vous n'avez pas). Vous avez juste à mettre la longitude/latitude de données ci-dessous le code dans un fichier texte appelé 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]-O-v[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))

Merci d'avance pour toutes suggestions ou des connaissances.

EDIT: Voici un schéma qui montre une projection de la sphère unité avec un polygone et le centre de gravité-je calculer à partir du code. Clairement, le centre de gravité est à tort que le polygone est plutôt petite et convexe, mais encore le centre de gravité se situe en dehors de son périmètre.

EDIT: Voici une très similaires ensemble de coordonnées à celles ci-dessus, mais dans l'original [lon,lat] format que j'ai l'habitude d'utiliser (qui est maintenant converti à [x,y,z] par le code mis à jour).

  -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

EDIT: quelques exemples...avec 4 sommets de la définition d'un carré parfait, centré à [1,0,0] - je obtenir le résultat attendu: Cependant, à partir d'un non-symétrique du triangle-je obtenir un centre de gravité qui n'est nulle part près...le centre de gravité tombe sur le côté le plus éloigné de la sphère (ici projeté sur le devant côté, comme l'antipode): Fait intéressant, le centre de gravité de l'estimation apparaît "stable" dans le sens que si je inverser la liste (à l'aller en sens horaire, antihoraire commande ou vice-versa), le centre de gravité est proportionnellement inverse exactement.

Answer 1

Je pense que cela va le faire. Vous devriez être en mesure de reproduire ce résultat il suffit de copier-coller le code ci-dessous.

• Vous aurez besoin d'avoir la latitude et la longitude de données dans un fichier appelé" longitude and latitude.txt. Vous pouvez copier-coller de l'original de l'échantillon de données est inclus ci-dessous le code.
• Si vous avez mplotlib il devra en outre produire le graphique ci-dessous
• Pour les non-évident calculs, j'ai inclus un lien qui explique ce qui se passe
• Dans le graphique ci-dessous, le vecteur de référence est très court (r = 1/10), de sorte que le 3d-centroïdes sont plus faciles à voir. Vous pouvez facilement retirer la mise à l'échelle afin de maximiser la précision.
• Note à l'op: j'ai réécrit presque tout donc je ne sais pas exactement où le code d'origine n'a pas de travail. Cependant, au moins je pense que c'était de ne pas prendre en considération la nécessité pour la poignée dans le sens horaire / antihoraire triangle de sommets.

---

Légende:

• (ligne noire) référence vecteur
• (petits points rouges) triangle sphérique 3d-centroïdes
• (gros rouge / bleu / vert dot) 3d-centre de gravité / projetées à la surface projetée sur le plan xy
• (bleu / vert lignes) de la partie sphérique de polygone et de la projection sur le plan 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 a-b-c 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 a-b-c
    # 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()

---

Voici la longitude et la latitude de données, vous pouvez les coller en 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

Je pense qu'une bonne approximation serait de calculer le centre de masse à l'aide d'pondérée en coordonnées cartésiennes et en projetant le résultat sur la sphère (en supposant que l'origine des coordonnées (0, 0, 0)^T).

Soit (p[0], p[1], ..., p[n-1]) les n points du polygone. L'env. (cartésien) centre de gravité peut être calculée par:

c = 1 / w * (somme de w[i] * p[i])

alors que w est la somme de tous les poids et tandis que le p[i] est un polygone point et w[i] est un poids pour ce point, par exemple

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

alors que |x| est la longueur d'un vecteur x. I. e. un point est pondérée avec la moitié de la longueur à la précédente et à la moitié de la longueur de la prochaine polygone point.

Ce centre de gravité c peut maintenant projetée sur la sphère par:

c' = r * c / |c|

alors que r est le rayon de la sphère.

Pour prendre en considération l'orientation du polygone (ccw cw), le résultat peut être

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