Mapping A Sphere To A Cube
This answer contains the cube2sphere and sphere2cube without the restriction of a = 1. So the cube has side 2a from -a to a and the radius of the sphere is a.
I know it's been 10 years since this question was asked. Nevertheless, I am giving the answer in case someone needs it. The implementation is in Python,
I am using (x, y, z) for the cube coordinates, (p, q, r) for the sphere coordinates and the relevant underscore variables (x_, y_, z_) meaning they have been produced by using the inverse function.
import math
from random import randint # for testing
def sign_aux(x):
return lambda y: math.copysign(x, y)
sign = sign_aux(1) # no built-in sign function in python, I know...
def cube2sphere(x, y, z):
if (all([x == 0, y == 0, z == 0])):
return 0, 0, 0
def aux(x, y_2, z_2, a, a_2):
return x * math.sqrt(a_2 - y_2/2 - z_2/2 + y_2*z_2/(3*a_2))/a
x_2 = x*x
y_2 = y*y
z_2 = z*z
a = max(abs(x), abs(y), abs(z))
a_2 = a*a
return aux(x, y_2, z_2, a, a_2), aux(y, x_2, z_2, a, a_2), aux(z, x_2, y_2, a, a_2)
def sphere2cube(p, q, r):
if (all([p == 0, q == 0, r == 0])):
return 0, 0, 0
def aux(s, t, radius):
A = 3*radius*radius
R = 2*(s*s - t*t)
S = math.sqrt( max(0, (A+R)*(A+R) - 8*A*s*s) ) # use max 0 for accuraccy error
iot = math.sqrt(2)/2
s_ = sign(s) * iot * math.sqrt(max(0, A + R - S)) # use max 0 for accuraccy error
t_ = sign(t) * iot * math.sqrt(max(0, A - R - S)) # use max 0 for accuraccy error
return s_, t_
norm_p, norm_q, norm_r = abs(p), abs(q), abs(r)
norm_max = max(norm_p, norm_q, norm_r)
radius = math.sqrt(p*p + q*q + r*r)
if (norm_max == norm_p):
y, z = aux(q, r, radius)
x = sign(p) * radius
return x, y, z
if (norm_max == norm_q):
z, x = aux(r, p, radius)
y = sign(q) * radius
return x, y, z
x, y = aux(p, q, radius)
z = sign(r) * radius
return x, y, z
# measuring accuracy
max_mse = 0
for i in range(100000):
x = randint(-20, 20)
y = randint(-20, 20)
z = randint(-20, 20)
p, q, r = cube2sphere(x, y, z)
x_, y_, z_ = sphere2cube(p, q, r)
max_mse = max(max_mse, math.sqrt(((x-x_)**2 + (y-y_)**2 + (z-z_)**2))/3)
print(max_mse)
# 1.1239159602905078e-07
max_mse = 0
for i in range(100000):
p = randint(-20, 20)
q = randint(-20, 20)
r = randint(-20, 20)
x, y, z = sphere2cube(p, q, r)
p_, q_, r_ = cube2sphere(x, y, z)
max_mse = max(max_mse, math.sqrt(((p-p_)**2 + (q-q_)**2 + (r-r_)**2))/3)
print(max_mse)
# 9.832883321715792e-08
Also, I mapped some points to check the function visually and these are the results.
After some rearranging you can get the "nice" forms
(1) 1/2 z^2 = (alpha) / ( y^2 - x^2) + 1
(2) 1/2 y^2 = (beta) / ( z^2 - x^2) + 1
(3) 1/2 x^2 = (gamma) / ( y^2 - z^2) + 1
where alpha = sx^2-sy^2 , beta = sx^2 - sz^2 and gamma = sz^2 - sy^2. Verify this yourself.
Now I neither have the motivation nor the time but from this point on its pretty straightforward to solve:
Substitute (1) into (2). Rearrange (2) until you get a polynomial (root) equation of the form
(4) a(x) * y^4 + b(x) * y^2 + c(x) = 0this can be solved using the quadratic formula for
y^2. Note thata(x),b(x),c(x)are some functions ofx. The quadratic formula yields 2 roots for (4) which you will have to keep in mind.Using (1),(2),(4) figure out an expression for
z^2in terms of onlyx^2.Using (3) write out a polynomial root equation of the form:
(5) a * x^4 + b * x^2 + c = 0where
a,b,care not functions but constants. Solve this using the quadratic formula. In total you will have 2*2=4 possible solutions forx^2,y^2,z^2pair meaning you will have 4*2=8 total solutions for possiblex,y,zpairs satisfying these equations. Check conditions on eachx,y,zpair and (hopefully) eliminate all but one (otherwise an inverse mapping does not exist.)
Good luck.
PS. It very well may be that the inverse mapping does not exist, think about the geometry: the sphere has surface area 4*pi*r^2 while the cube has surface area 6*d^2=6*(2r)^2=24r^2 so intuitively you have many more points on the cube that get mapped to the sphere. This means a many to one mapping, and any such mapping is not injective and hence is not bijective (i.e. the mapping has no inverse.) Sorry but I think you are out of luck.
----- edit --------------
if you follow the advice from MO, setting z=1 means you are looking at the solid square in the plane z=1.
Use your first two equations to solve for x,y, wolfram alpha gives the result:
x = (sqrt(6) s^2 sqrt(1/2 (sqrt((2 s^2-2 t^2-3)^2-24 t^2)+2 s^2-2 t^2-3)+3)-sqrt(6) t^2 sqrt(1/2 (sqrt((2 s^2-2 t^2-3)^2-24 t^2)+2 s^2-2 t^2-3)+3)-sqrt(3/2) sqrt((2 s^2-2 t^2-3)^2-24 t^2) sqrt(1/2 (sqrt((2 s^2-2 t^2-3)^2-24 t^2)+2 s^2-2 t^2-3)+3)+3 sqrt(3/2) sqrt(1/2 (sqrt((2 s^2-2 t^2-3)^2-24 t^2)+2 s^2-2 t^2-3)+3))/(6 s)
and
y = sqrt(-sqrt((2 s^2-2 t^2-3)^2-24 t^2)-2 s^2+2 t^2+3)/sqrt(2)
where above I use s=sx and t=sy, and I will use u=sz. Then you can use the third equation you have for u=sz. That is lets say that you want to map the top part of the sphere to the cube. Then for any 0 <= s,t <= 1 (where s,t are in the sphere's coordinate frame ) then the tuple (s,t,u) maps to (x,y,1) (here x,y are in the cubes coordinate frame.) The only thing left is for you to figure out what u is. You can figure this out by using s,t to solve for x,y then using x,y to solve for u.
Note that this will only map the top part of the cube to only the top plane of the cube z=1. You will have to do this for all 6 sides (x=1, y=1, z=0 ... etc ). I suggest using wolfram alpha to solve the resulting equations you get for each sub-case, because they will be as ugly or uglier as those above.
I want to give gmatt credit for this because he's done a lot of the work. The only difference in our answers is the equation for x.
To do the inverse mapping from sphere to cube first determine the cube face the sphere point projects to. This step is simple - just find the component of the sphere vector with the greatest length like so:
// map the given unit sphere position to a unit cube position
void cubizePoint(Vector3& position) {
double x,y,z;
x = position.x;
y = position.y;
z = position.z;
double fx, fy, fz;
fx = fabsf(x);
fy = fabsf(y);
fz = fabsf(z);
if (fy >= fx && fy >= fz) {
if (y > 0) {
// top face
position.y = 1.0;
}
else {
// bottom face
position.y = -1.0;
}
}
else if (fx >= fy && fx >= fz) {
if (x > 0) {
// right face
position.x = 1.0;
}
else {
// left face
position.x = -1.0;
}
}
else {
if (z > 0) {
// front face
position.z = 1.0;
}
else {
// back face
position.z = -1.0;
}
}
}
For each face - take the remaining cube vector components denoted as s and t and solve for them using these equations, which are based on the remaining sphere vector components denoted as a and b:
s = sqrt(-sqrt((2 a^2-2 b^2-3)^2-24 a^2)+2 a^2-2 b^2+3)/sqrt(2)
t = sqrt(-sqrt((2 a^2-2 b^2-3)^2-24 a^2)-2 a^2+2 b^2+3)/sqrt(2)
You should see that the inner square root is used in both equations so only do that part once.
Here's the final function with the equations thrown in and checks for 0.0 and -0.0 and the code to properly set the sign of the cube component - it should be equal to the sign of the sphere component.
void cubizePoint2(Vector3& position)
{
double x,y,z;
x = position.x;
y = position.y;
z = position.z;
double fx, fy, fz;
fx = fabsf(x);
fy = fabsf(y);
fz = fabsf(z);
const double inverseSqrt2 = 0.70710676908493042;
if (fy >= fx && fy >= fz) {
double a2 = x * x * 2.0;
double b2 = z * z * 2.0;
double inner = -a2 + b2 -3;
double innersqrt = -sqrtf((inner * inner) - 12.0 * a2);
if(x == 0.0 || x == -0.0) {
position.x = 0.0;
}
else {
position.x = sqrtf(innersqrt + a2 - b2 + 3.0) * inverseSqrt2;
}
if(z == 0.0 || z == -0.0) {
position.z = 0.0;
}
else {
position.z = sqrtf(innersqrt - a2 + b2 + 3.0) * inverseSqrt2;
}
if(position.x > 1.0) position.x = 1.0;
if(position.z > 1.0) position.z = 1.0;
if(x < 0) position.x = -position.x;
if(z < 0) position.z = -position.z;
if (y > 0) {
// top face
position.y = 1.0;
}
else {
// bottom face
position.y = -1.0;
}
}
else if (fx >= fy && fx >= fz) {
double a2 = y * y * 2.0;
double b2 = z * z * 2.0;
double inner = -a2 + b2 -3;
double innersqrt = -sqrtf((inner * inner) - 12.0 * a2);
if(y == 0.0 || y == -0.0) {
position.y = 0.0;
}
else {
position.y = sqrtf(innersqrt + a2 - b2 + 3.0) * inverseSqrt2;
}
if(z == 0.0 || z == -0.0) {
position.z = 0.0;
}
else {
position.z = sqrtf(innersqrt - a2 + b2 + 3.0) * inverseSqrt2;
}
if(position.y > 1.0) position.y = 1.0;
if(position.z > 1.0) position.z = 1.0;
if(y < 0) position.y = -position.y;
if(z < 0) position.z = -position.z;
if (x > 0) {
// right face
position.x = 1.0;
}
else {
// left face
position.x = -1.0;
}
}
else {
double a2 = x * x * 2.0;
double b2 = y * y * 2.0;
double inner = -a2 + b2 -3;
double innersqrt = -sqrtf((inner * inner) - 12.0 * a2);
if(x == 0.0 || x == -0.0) {
position.x = 0.0;
}
else {
position.x = sqrtf(innersqrt + a2 - b2 + 3.0) * inverseSqrt2;
}
if(y == 0.0 || y == -0.0) {
position.y = 0.0;
}
else {
position.y = sqrtf(innersqrt - a2 + b2 + 3.0) * inverseSqrt2;
}
if(position.x > 1.0) position.x = 1.0;
if(position.y > 1.0) position.y = 1.0;
if(x < 0) position.x = -position.x;
if(y < 0) position.y = -position.y;
if (z > 0) {
// front face
position.z = 1.0;
}
else {
// back face
position.z = -1.0;
}
}
So, this solution isn't nearly as pretty as the cube to sphere mapping, but it gets the job done!
Any suggestions to improve the efficiency or read ability of the code above are appreciated!
--- edit --- I should mention that I have tested this and so far in my tests the code appears correct with the results being accurate to at least the 7th decimal place. And that was from when I was using floats, it's probably more accurate now with doubles.
--- edit --- Here's an optimized glsl fragment shader version by Daniel to show that it doesn't have to be such a big scary function. Daniel uses this to filter sampling on cube maps! Great idea!
const float isqrt2 = 0.70710676908493042;
vec3 cubify(const in vec3 s)
{
float xx2 = s.x * s.x * 2.0;
float yy2 = s.y * s.y * 2.0;
vec2 v = vec2(xx2 – yy2, yy2 – xx2);
float ii = v.y – 3.0;
ii *= ii;
float isqrt = -sqrt(ii – 12.0 * xx2) + 3.0;
v = sqrt(v + isqrt);
v *= isqrt2;
return sign(s) * vec3(v, 1.0);
}
vec3 sphere2cube(const in vec3 sphere)
{
vec3 f = abs(sphere);
bool a = f.y >= f.x && f.y >= f.z;
bool b = f.x >= f.z;
return a ? cubify(sphere.xzy).xzy : b ? cubify(sphere.yzx).zxy : cubify(sphere);
}