Better GPU procedural generated Spherified Cube Sphere

better_spherified_cube_sphere.glsl

// Code (C) by Benjamin 'BeRo' Rosseaux - zlib licensed

// Based on the ideas from Stephen Cameron's https://scaryreasoner.wordpress.com/2016/01/23/thoughts-on-tesselating-a-sphere/ blog post.

// ...

void getCubeSphereNormals(const in int face, in vec4 uv0011, out vec3 vectors[4]){
  const float deltaAngle = 1.5707963267948966; // radians(90.0);
  const float startAngle = 0.7853981633974483; // radians(45.0);
  const vec3 normals[6] = vec3[6](
    vec3(1.0, 0.0, 0.0), // +X
    vec3(-1.0, 0.0, 0.0), // -X
    vec3(0.0, 1.0, 0.0), // +Y
    vec3(0.0, -1.0, 0.0), // -Y
    vec3(0.0, 0.0, 1.0), // +Z
    vec3(0.0, 0.0, -1.0) // -Z
  );
  const ivec2 angleMap[6] = ivec2[6]( // x = angle A axis, y = angle B axis
    ivec2(2, 1), // +X
    ivec2(2, 1), // -X
    ivec2(0, 2), // +Y
    ivec2(0, 2), // -Y
    ivec2(0, 1), // +Z
    ivec2(0, 1)  // -Z
  );
  const vec4 signMap[6] = vec4[6]( // x = start angle A sign, y = start angle B sign, z = delta angle A sign, w = delta angle B sign
    vec4(1.0, 1.0, -1.0, -1.0), // +X
    vec4(1.0, -1.0, -1.0, 1.0), // -X
    vec4(1.0, 1.0, -1.0, -1.0), // +Y
    vec4(1.0, -1.0, -1.0, 1.0), // -Y
    vec4(1.0, -1.0, -1.0, 1.0), // +Z
    vec4(-1.0, -1.0, 1.0, 1.0)  // -Z
  );   
  const vec4 angles = tan(fma(fma(uv0011, vec2(-1.0, 1.0).xyxy, vec2(1.0, 0.0).xyxy) * deltaAngle, signMap[face].zwzw, signMap[face].xyxy * startAngle));
  const ivec2 angleIndices = angleMap[face];
  const vec3 baseBormal = normals[face];
  [[unroll]] for(uint quadVertexIndex = 0u; quadVertexIndex < 4u; quadVertexIndex++){
    vectors[quadVertexIndex] = baseBormal;
    vectors[quadVertexIndex][angleIndices.x] = angles[uvec4(0u, 2u, 2u, 0u)[quadVertexIndex]];
    vectors[quadVertexIndex][angleIndices.y] = angles[uvec4(1u, 1u, 3u, 3u)[quadVertexIndex]];
    vectors[quadVertexIndex] = normalize(vectors[quadVertexIndex]);
  }
}

// ...

    uint vertexIndex = uint(gl_VertexIndex);

    // 0xe24 = 3,2,0,2,1,0 (two bit wise encoded triangle indices, reversed for bitshifting for 0,1,2, 0,2,3 output order)
    // 0xb4 = 180 = 0b10110100 (bitwise encoded x y coordinates, where x is the first bit and y is the second bit in every two-bit pair)   

    uint quadIndex = vertexIndex / 6u,   
         quadVertexIndex = (0xe24u >> ((vertexIndex - (quadIndex * 6u)) << 1u)) & 3u, 
         quadVertexUVIndex = (0xb4u >> (quadVertexIndex << 1u)) & 3u;     

    vec3 vectors[4];
    getCubeSphereNormals(int(sideIndex), vec4(uvec2(sideQuadX, sideQuadY).xyxy + uvec4(0u, 0u, 1u, 1u)) / vec4(resolution), vectors);

    // Find the two triangles by the shortest diagonal and adjust quadVertexUVIndex accordingly.
    if(distance(vectors[1], vectors[3]) < distance(vectors[0], vectors[2])){
      quadVertexUVIndex = uvec4(1u, 3u, 0u, 2u)[quadVertexUVIndex];
    }

    sphereNormal = vectors[uvec4(0u, 1u, 3u, 2u)[quadVertexUVIndex]];

    // Fast cheap orthonormal basis as tangent space
    vec3 sphereTangent = normalize(cross((abs(sphereNormal.y) < 0.999999) ? vec3(0.0, 1.0, 0.0) : vec3(0.0, 0.0, 1.0), sphereNormal));
    vec3 sphereBitangent = normalize(cross(sphereNormal, sphereTangent));

// ...