ohhh now I see what you mean.
You have a function going from the position on the unit cube C to the unit sphere S and you would like a function going the other way round. I got it that far, what I didn't understand is that, since the parameterization of the sphere you use is not the standard "gnomonic" one (the standard the cubemap), you cannot use the standard inverse function, which both me and Polly discussed.
I can see three practical solutions.
Way 1) Go back to the default "vanilla" cubemap.
I.e. keep things simple. It is true that the parameterization you use is more isometric (preserves better both angles and areas, meaning the squares are more regularly sized and shaped on the sphere), but ask yourself, is it worth it. A terrain would probably look very similar.
Way 2) Only ever use cube->sphere (never sphere->cube directly).
Like this: you have a point p = (x,y,z) on the cube and with your function you can find s = (sx,sy,sz) on the sphere. If you want to find neighboring points, move p over the cube, get p'=(x',y',z'), and project it on the sphere, funding s' = (sx',sy',sz') which will be close enough to s.
Way 3) Do the math and compute sphere->cube
I think you'll like this one more, but it takes a little more effort...
I did the math for you (initial, and not tested, but it should be close).
So you have a point (sx, sy, sz) on the sphere and you want to find
first, the face of the cube (one of -X, +X, -Y, +Y, -Z, +Z) and then a position on that face.
So basically you want to find (x,y,z) on the cube, i.e. such that one of x,y,z is -1 or +1, and the other two coords are between -1 and +1 (included). Ok?
Ok, first of all: in your math above, your sphere seems to have unit DIAMETER (i.e. 0.5 radius). This is a bit of an useless complication. It is easier if you remove that "/2" at the end of sx = ... (this is valid both in your parametrization and in standard parametrization). Once you have your position on the radius-one sphere, just multiply by planet radius.
So your formulas look like this:
sx = Sqrt(1.0 - yy/2 - zz/2 + ( (yy * zz) / 3.0 ) ) * x (1)
sy = Sqrt(1.0 - zz/2 - xx/2 + ( (xx * zz) / 3.0 ) ) * y (2)
sz = Sqrt(1.0 - xx/2 - yy/2 + ( (xx * yy) / 3.0 ) ) * z (3)
(I numbered the three formulas)
Where
(sx,sy,sz) = position on the UNIT sphere (zero centered, diameter two)
(x,y,z) = position on the UNIT cube (zero centered, side two)
and naturally xx = x*x
Now, given (sx,sy,sz) you'll find (x,y,z) on the cube.
I.e. which face and position on that face.
It is not difficult.
Step 1: determine the face. I.e. determine which of x,y,z is +1 or -1.
This is the same as before. Grab the biggest between |sx|, |sy|, |sz|, and see the sign before the abs:
For example, if |sx| is then biggest then: if sx<0 then x=-1, if sx>0 then x=+1.
Ok? Probably your code will have 6 cases, x=+1, x=-1, y=+1... etc.
Say it was x=+1 (i.e. is the +X face).
The other cases will be similar.
Next, you need to determine y and z (which will be the position inside that face).
Here is how to do so.
Just take (2) an (3) (you don't care for 1, because you know x already) and substitute x=+1:
you get
sy = Sqrt(0.5 - zz/6 ) * y
sz = Sqrt(0.5 - yy/6 ) * z
(this will not be your code just yet, it is just a way to get there)
That's a lot simpler already, right?
Now let's invert it. You get (after a few passages):
syy = (0.5 - zz/6) * yy
szz = (0.5 - yy/6) * zz
where syy = sy*sy and szz = sz*sz.
Going further, we get
yy = syy/(0.5-zz/6) = 6*syy/(3-zz)
and
0.5*zz*zz + zz*(syy-szz-1.5) + 3*szz = 0
That is degree 4 in z (zz*zz is z^4) but fear not.
We will just solve it for zz (i.e. our unkown is zz, not z)
Now, just rename :
a = 0.5
b = (syy-szz-1.5)
c = 3*szz
that becomes
a*zz*zz + b*zz + c = 0
Which is solved by zz = -b +sqrt( b*b-2*c). // because a = 0.5
(you'll need the +sqrt solution, although the solution with -sqrt is also a solution of the equation).
Now you have zz and you can find yy, as
yy = 6*syy/(3-zz)
Last passage. Now you have yy and zz, the squares of the z and y you are looking for.
So, should you take y = sqrt(yy) or y = -sqrt(yy)?
Luckily, that is easy to determine: if you think about it,
the sign of y is the same as the sign of sy (and same for z).
So in total the final computation for this case is:
function find_pos_on_cube( float sx, sy, sz )
{
// case |sx| is biggest, and sx>0.
// Face is: +X.
// Position on the cube: x = +1.
//
float syy = sy*sy;
float szz = sz*sz;
float b = (syy-szz-1.5);
float c = 3*szz;
float zz = -b + sqrt( b*b - 2*c);
float yy = 6*syy/(3-zz);
float x,y,z; // pos on cube
x = +1;
z = (sz>0)? sqrt(zz) : -sqrt(zz);
y = (sy>0)? sqrt(yy) : -sqrt(yy);
}
Developer
