Issue 1:
When all of the points are in one hemisphere some of the Voronoi points are incorrect.

To be more precise, certain of the points are on the wrong side of the sphere. An example is to create a tetrahedron with the following generator points:

generators (polar)

[[-25.56198973 10.81015471]

[-54.40661695 26.73672285]

[-55.29180524 -6.7928405 ]

[-45.12839961 9.12243389]]

generators (rectangular)

[[ 0.88610999 -0.42383021 0.18755541]

[ 0.51980039 -0.72622668 0.4498915 ]

[ 0.56540011 -0.81629197 -0.11827989]

[ 0.69659682 -0.69972598 0.15854467]]

The first three Voronoi points and associated Delaunay triangles come out fine:

Voronoi vertex number 1 , coordinates: [ 0.73647759 -0.54764287 0.39709955]

Delaunay vertices

[[ 0.69659682 -0.69972598 0.15854467]

[ 0.88610999 -0.42383021 0.18755541]

[ 0.51980039 -0.72622668 0.4498915 ]]

Voronoi vertex number 2 , coordinates: [ 0.42155897 -0.88977832 0.17487877]

Delaunay vertices

[[ 0.69659682 -0.69972598 0.15854467]

[ 0.56540011 -0.81629197 -0.11827989]

[ 0.51980039 -0.72622668 0.4498915 ]]

Voronoi vertex number 3 , coordinates: [ 0.82153008 -0.54761634 -0.15875979]

Delaunay vertices

[[ 0.69659682 -0.69972598 0.15854467]

[ 0.56540011 -0.81629197 -0.11827989]

[ 0.88610999 -0.42383021 0.18755541]]

The fourth point is not correct - it should be on the opposite side of the sphere:

Voronoi vertex number 0 , coordinates: [ 0.69593711 -0.69851441 0.16658079]

Delaunay vertices

[[ 0.56540011 -0.81629197 -0.11827989]

[ 0.88610999 -0.42383021 0.18755541]

[ 0.51980039 -0.72622668 0.4498915 ]]

I have a good understanding of the algorithm and what the problem is but I don't know how to fix it using qhull.

Issue 2:
As pointed out by @tylerjereddy in #5232, the scipy.spatial.SphericalVoronoi algorithm does not cover cases where all generator points are on a single hemisphere. How would an extension of the algorithm look like to include those cases? Are there theoretical treatises tackling these cases?

Reproducing code example:

import scipy.spatial
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import proj3d
from matplotlib import colors
from mpl_toolkits.mplot3d.art3d import Poly3DCollection

# calculate points in single hemisphere
phi = np.linspace(0,2*np.pi,10,endpoint=False) # azimuth angle
theta = np.linspace(0.001,np.pi*0.4,5) # polar angle
theta = theta[np.newaxis, :].T
radius = 1;

phiv, thetav = np.meshgrid(phi, theta)
phiv = np.reshape(phiv,(50,1))
thetav = np.reshape(thetav,(50,1))

x = radius*np.cos(phiv)*np.sin(thetav)
y = radius*np.sin(phiv)*np.sin(thetav)
z = radius*np.cos(thetav)

points = np.concatenate([x,y,z],axis=1)

# Voronoi-Regions
sv = scipy.spatial.SphericalVoronoi(points)
sv.sort_vertices_of_regions()

# plotting
fig = plt.figure()
ax = fig.add_subplot('111', projection='3d')
ax.scatter(points[...,0], points[...,1], points[...,2])

for region in sv.regions:
random_color = colors.rgb2hex(np.random.rand(3))
polygon = Poly3DCollection([sv.vertices[region]], alpha=1.0)
polygon.set_color(random_color)
ax.add_collection3d(polygon)

ax.set_xlim([-1.0,1.0])
ax.set_ylim([-1.0,1.0])
ax.set_zlim([-1.0,1.0])
ax.view_init(45,0)

fig.savefig('hemisphere.png', dpi=300)

resulting plot

Scipy/Numpy/Python version information:

1.1.0 1.14.3 sys.version_info(major=3, minor=5, micro=2, releaselevel='final', serial=0)