convex_hull does not check the points are in one hemisphere

[anutkk]

The function convex_hull in polygon.py does not check that all points are in one open hemisphere. If the points cover an area greater than an open hemisphere, the convex hull should be the entire sphere. Instead, the function returns part of the points. The solution is to check for "hemisphericity" before applying the modified Graham's algorithm.

The code:

import os
import spherical_geometry
from spherical_geometry import polygon
import random
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from itertools import product, combinations

lons = np.array([])
lats = np.array([])
number_of_points = 20

for jj in range(0, number_of_points):
    u = random.random()
    v = random.random()
    # convert to spherical coordinates
    theta = 2 * np.pi * u * 0.5 + np.pi # range [0,2pi)
    phi = 2*np.arccos(2 * v - 1) # range (0, pi/2]
    # convert to lng,lat
    lng = theta / (2 * np.pi) * 360 - 180
    lat =   phi / (2 * np.pi) * 360 - 90
    lons = np.append(lons, lng)
    lats = np.append(lats, lat)
R=1
lonr = np.deg2rad(lons)
latr = np.deg2rad(lats)
pts_x = R * np.cos(latr)*np.cos(lonr)
pts_y = R * np.cos(latr)*np.sin(lonr)
pts_z = R * np.sin(latr)
a = np.column_stack((pts_x,pts_y,pts_z))
aa = spherical_geometry.polygon.SphericalPolygon.convex_hull(a)
bb=aa.polygons[0]
cc=bb.points

fig = plt.figure()
ax = fig.gca(projection='3d')

earthRadius =1 # 6378.137;

# draw sphere
u, v = np.mgrid[0:2*np.pi:20j, 0:np.pi:10j]
x = earthRadius * np.cos(u)*np.sin(v)
y = earthRadius * np.sin(u)*np.sin(v)
z = earthRadius * np.cos(v)
ax.plot_wireframe(x, y, z, color="k")


ax.scatter(pts_x, pts_y, pts_z, color="g", s=50)
ax.plot(cc[:,0], cc[:,1], cc[:,2], color="r")

plt.show()

Disclosure: I maintain a (new and for now lacking a lot of functions) parallel package of spherical geometry with a different focus.