Mapping a Cube to a Sphere

I have been looking for a way to map a unit cube (with vertices $x^{2}=1$, $y^{2}=1$, $z^{2}=1$) to a unit sphere ($x^{2}+y^{2}+z^{2}=1$) with minimal distortion of the great circles formed by mapping the coordinate lines on the cube face. As can be seen from the following picture, a simple radial contraction onto the sphere surface leads to large visible distortion of the great circles. 

Phil Nowell (mapping-cube-to-sphere) derived an elegant mapping that generates a much more uniform subdivision of the great circles.

However, close inspection shows that there is still some room for improvement. Points near the center of the cube faces get more compressed than those near the edges.

By using the rotation of central planes to map points on the cube to equidistant points on the sphere, I was able to come up with the following expressions for the coordinates of the sphere:

$x_{sphere}=\frac{x_{c}}{\sqrt{x^{2}_{c}+y^{2}_{c}+z^{2}_{c}} } $

$y_{sphere}=\frac{y_{c}}{\sqrt{x^{2}_{c}+y^{2}_{c}+z^{2}_{c}} } $

$z_{sphere}=\frac{z_{c}}{\sqrt{x^{2}_{c}+y^{2}_{c}+z^{2}_{c}} } $

where

$x_{c}=\sqrt{\left( x^{p}+y^{2}+z^{2}\right)  } tan\left( x\  arctan\left( \frac{1}{\sqrt{\left( x^{p}+y^{2}+z^{2}\right)  } } \right)  \right) $

$y_{c}=\sqrt{\left( x^{2}+y^{p}+z^{2}\right)  } tan\left( y\  arctan\left( \frac{1}{\sqrt{\left( x^{2}+y^{p}+z^{2}\right)  } } \right)  \right) $

$z_{c}=\sqrt{\left( x^{2}+y^{2}+z^{p}\right)  } tan\left( z\  arctan\left( \frac{1}{\sqrt{\left( x^{2}+y^{2}+z^{p}\right)  } } \right)  \right) $

p is a large even number and

x, y, z are the coordinates of the cube.

The picture below is obtained with p=50

To express the quality of the mapping, I calculate the ratio of the maximum and minimum length of a line segment on the most distorted great circle, i, as follows:

$$Q=\frac{l_{i,min}}{l_{i,max}} \cdot 100\% $$

so that a perfect distribution gives a quality of 100%

With this, the quality of the radial contraction of the 9x9x9 grid depicted above becomes 58%, that of Nowell’s mapping becomes 79% and my proposed mapping (with p=50) becomes 96%. It is also worth noting that the quality on the central coordinate lines and the edges is 100%.

Update 1

The following iterative scheme gives a further improvement:

$x_{c,0} = x$,  $y_{c,0} = y$,  $z_{c,0} = z$

$\begin{gathered}x_{c,i+1}=\sqrt{x^{p}_{c,i}+y^{2}_{c,i}+z^{2}_{c,i}} tan(x\  arctan(\frac{1}{\sqrt{x^{p}_{c,i}+y^{2}_{c,i}+z^{2}_{c,i}} } ))\\ y_{c,i+1}=\sqrt{x^{2}_{c,i}+y^{p}_{c,i}+z^{2}_{c,i}} tan(y\  arctan(\frac{1}{\sqrt{x^{2}_{c,i}+y^{p}_{c,i}+z^{2}_{c,i}} } ))\\ z_{c,i+1}=\sqrt{x^{2}_{c,i}+y^{2}_{c,i}+z^{p}_{c,i}} tan(z\  arctan(\frac{1}{\sqrt{x^{2}_{c,i}+y^{2}_{c,i}+z^{p}_{c,i}} } ))\end{gathered} $

With this, the quality of the mapped 9x9x9 grid increases after 2 iterations from 96.1% to 97.7%

For finer grids the improvement is more pronounced. A 50x50x50 grid (with p=500) achieves an improvement from 95.3% to 97.1% after 1 iteration and a grid of 500x500x500 (with p=5000) shows an improvement from 95.2% to 97.0% with 1 iteration

Background

With reference to the below figure, equidistant points on a great circle can be obtained by subdividing angle $\phi_{max}$ in equal parts. By using the coordinate on the cube as a scale factor for $\phi_{max}$, equidistant points on the cube will  be mapped to equidistant points on the great circle.

 

The process of mapping the original point on the cube to the sphere can thus be broken down in two steps, i.e. a mapping of the original point on the cube to an image point on the cube and the contraction of the latter point to the sphere. The coordinates of the image point on the cube ($x_{c}, y_{c}, z_{c}$) can be obtained by extending the radius of the image on the sphere out to the surface of the cube. In this way the following equation is derived:

$x_{c}=\sqrt{\left(y^{2}+z^{2}\right)  } tan\left( x\  arctan\left( \frac{1}{\sqrt{\left( y^{2}+z^{2}\right)  } } \right)  \right) $

To avoid a singularity at y=z=0, while avoiding a test for this singularity, I add the term $x^p$, which, for large even p, tends to 0 when $\left\vert \text{} x\right\vert <1$ and to 1 for $\left\vert \text{} x\right\vert=1$ (which is the case on the cube at y=z=0). Having said that, this singularity can also be handled in other ways.