Cube to Sphere mapping (inverse function wanted)

Question

I've came across a cube to sphere mapping function that provides a more uniform result than just normalizing the coordinates or other mapping methods. Unfortunately there is no unwrapping function.

Source: http://mathproofs.blogspot.com/2005/07/mapping-cube-to-sphere.html

vec3 spherify ( vec3 v ) {
    float x2 = v.x * v.x;
    float y2 = v.y * v.y;
    float z2 = v.z * v.z;
    vec3 s;
    s.x = v.x * sqrt(1.0 - y2 / 2.0 - z2 / 2.0 + y2 * z2 / 3.0);
    s.y = v.y * sqrt(1.0 - x2 / 2.0 - z2 / 2.0 + x2 * z2 / 3.0);
    s.z = v.z * sqrt(1.0 - x2 / 2.0 - y2 / 2.0 + x2 * y2 / 3.0);
    return s;
}

How could this be unwrapped back to a cube face? To wrap back to a square i use the following, though it doesn't unwrap the special squeezing of the coordinates.

vec3 cubify ( vec3 s ) {
    
    s.x = ( s.x / s.z );
    s.y = ( s.y / s.z );
    
    return s;
}

Which looks as the following Mapping

Answer 1

How about something like this (C++/VCL):

//---------------------------------------------------------------------------
#include <vcl.h>
#pragma hdrstop
#include "Unit1.h"
#include "gl_simple.h"
#include "glsl_math.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
const int n=10*10*6;
vec3 col[n];
vec3 cube[n];
vec3 sphere[n];
vec3 cube2[n];
//---------------------------------------------------------------------------
vec3 spherify(vec3 v)
    {
    float x2 = v.x * v.x;
    float y2 = v.y * v.y;
    float z2 = v.z * v.z;
    vec3 s;
    s.x = v.x * sqrt(1.0 - y2 / 2.0 - z2 / 2.0 + y2 * z2 / 3.0);
    s.y = v.y * sqrt(1.0 - x2 / 2.0 - z2 / 2.0 + x2 * z2 / 3.0);
    s.z = v.z * sqrt(1.0 - x2 / 2.0 - y2 / 2.0 + x2 * y2 / 3.0);
    return s;
    }
//---------------------------------------------------------------------------
vec3 cubify(vec3 v)
    {
    int i;
    float r,a;
    // major axis and size
    a=fabs(v.x);         { r=a; i=0; }
    a=fabs(v.y); if (r<a){ r=a; i=1; }
    a=fabs(v.z); if (r<a){ r=a; i=2; }
    v/=r; r*=1.75; a=4.0*r/M_PI;
    // convert of cube + linearization
         if (i==0){ v.y=a*atan(v.y/r); v.z=a*atan(v.z/r); }
    else if (i==1){ v.x=a*atan(v.x/r); v.z=a*atan(v.z/r); }
    else          { v.x=a*atan(v.x/r); v.y=a*atan(v.y/r); }
    // just remedy boundaries after linearization
    if (v.x<-1.0) v.x=-1.0;
    if (v.x>+1.0) v.x=+1.0;
    if (v.y<-1.0) v.y=-1.0;
    if (v.y>+1.0) v.y=+1.0;
    if (v.z<-1.0) v.z=-1.0;
    if (v.z>+1.0) v.z=+1.0;
    return v;
    }
//---------------------------------------------------------------------------
void set_cube()
    {
    float u,v,d;
    int m=sqrt(n/6),i,j,k;
    k=0; d=2.0/float(m-1);
    for (u=-1.0,i=0;i<m;i++,u+=d)
     for (v=-1.0,j=0;j<m;j++,v+=d)
        {
        col[k]=vec3(0.5,0.0,0.0); cube[k]=vec3(u,v,-1.0); k++;
        col[k]=vec3(1.0,0.0,0.0); cube[k]=vec3(u,v,+1.0); k++;
        col[k]=vec3(0.0,0.5,0.0); cube[k]=vec3(u,-1.0,v); k++;
        col[k]=vec3(0.0,1.0,0.0); cube[k]=vec3(u,+1.0,v); k++;
        col[k]=vec3(0.0,0.0,0.5); cube[k]=vec3(-1.0,u,v); k++;
        col[k]=vec3(0.0,0.0,1.0); cube[k]=vec3(+1.0,u,v); k++;
        }
    }
//---------------------------------------------------------------------------
void gl_draw()
    {
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    float aspect=float(xs)/float(ys);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    gluPerspective(60.0/aspect,aspect,0.1,100.0);
    glMatrixMode(GL_TEXTURE);
    glLoadIdentity();
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    glTranslatef(0.0,0.0,-10.5);
    glRotatef(-10.0,1.0,0.0,0.0);
    glRotatef(-20.0,0.0,1.0,0.0);

    glEnable(GL_DEPTH_TEST);
    glDisable(GL_TEXTURE_2D);

    int i;

    glPointSize(2);

    glMatrixMode(GL_MODELVIEW);
    glTranslatef(-3.0,0.0,0.0); glBegin(GL_POINTS); for (i=0;i<n;i++){ glColor3fv(col[i].dat); glVertex3fv(cube[i].dat); } glEnd();     // set_cube
    glTranslatef(+3.0,0.0,0.0); glBegin(GL_POINTS); for (i=0;i<n;i++){ glColor3fv(col[i].dat); glVertex3fv(sphere[i].dat); } glEnd();   // spherify
    glTranslatef(+3.0,0.0,0.0); glBegin(GL_POINTS); for (i=0;i<n;i++){ glColor3fv(col[i].dat); glVertex3fv(cube2[i].dat); } glEnd();    // cubify

    glEnd();
    glPointSize(1);

    glFlush();
    SwapBuffers(hdc);
    }
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner):TForm(Owner)
    {
    gl_init(Handle);
    int i;
    set_cube();
    for (i=0;i<n;i++) sphere[i]=spherify(cube[i]);
    for (i=0;i<n;i++)  cube2[i]=cubify(sphere[i]);
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormDestroy(TObject *Sender)
    {
    gl_exit();
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormPaint(TObject *Sender)
    {
    gl_draw();
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::Timer1Timer(TObject *Sender)
    {
    gl_draw();
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormResize(TObject *Sender)
    {
    gl_resize(ClientWidth,ClientHeight);
    gl_draw();
    }
//---------------------------------------------------------------------------

Just ignore the VCL stuff. Code creates uniform grid cube points using set_cube , that is converted into sphere using your spherify and that is finally converted to cube2 using mine cubify.

Here preview:

from left cube,sphere,cube2. The colors are stored in col to better show the mapping between points...

The idea behind cubify is to leave biggest coordinate as is and the other two convert into spherical angle and then use this angle as coordinate. Basically its a reverse of this. Its a bit nonlinear near edges hence the slight shifting of 45 deg range to smaller ones... Also to avoid crossings of points above surface after linearization another check is involved (the 6 ifs at the end).