Recursion-free procedural GPU icosphere

icosphere.glsl

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

// ...

#extension GL_EXT_control_flow_attributes : enable

// ...

uint resolution = pushConstants.resolution;
uint squaredResolution = pushConstants.resolution * pushConstants.resolution;
uint countTotalVertices = squaredResolution * 3u * 20u; // 3 vertices per triangle, 20 triangles per icosahedron

// ...

mat3 tessellateTriangle(const in uint index, const in uint resolution, const in vec3 inputVertices[3], out vec3 outputVertices[3]){

  // Setup some variables for the barycentric coordinates
  const float y = floor(sqrt(float(index))), 
              x = (float(index) - (y * y)) * 0.5, 
              inverseResolution = 1.0 / float(resolution);

  // Check if it is a inverted triangle case.
  const bool inversed = fract(x) > 0.4;
  
  // Calculate the barycentric coordinates of the triangle, which is made of three points.
  const vec2 bc0 = (vec2(x, float(resolution) - y) * inverseResolution) + (inversed ? vec2(0.5 * inverseResolution, 0.0) : vec2(0.0));
  const vec2 bc1 = (bc0.xy + vec2(inverseResolution, -inverseResolution)) - (inversed ? vec2(inverseResolution, 0.0) : vec2(0.0));
  const vec2 bc2 = bc0.xy - (inversed ? vec2(inverseResolution, 0.0) : vec2(0.0, inverseResolution));
  
  // Put the barycentric coordinates into a 3x3 matrix for easier access, including the third w coordinate, which is just 1.0 - (u + v).
  mat3 result = mat3(
    vec3(bc0.xy, 1.0 - (bc0.x + bc0.y)), 
    vec3(bc1.xy, 1.0 - (bc1.x + bc1.y)), 
    vec3(bc2.xy, 1.0 - (bc2.x + bc2.y))
  );

  // Maybe not really necessary, but just for safety reasons, clamp the barycentric coordinates to the triangle for to avoid possible out-of-bound coordinates.
  [[unroll]] for(uint barycentricIndex = 0u; barycentricIndex < 3u; barycentricIndex++){
    vec3 uvw = result[barycentricIndex];
    vec3 p = (inputVertices[0] * uvw.x) + (inputVertices[1] * uvw.y) + (inputVertices[2] * uvw.z);
    if(uvw.x < 0.0){
      float t = clamp(dot(p - inputVertices[1], inputVertices[2] - inputVertices[1]) / dot(inputVertices[2] - inputVertices[1], inputVertices[2] - inputVertices[1]), 0.0, 1.0);
      result[barycentricIndex] = vec3(0.0, 1.0 - t, t);
    }else if(uvw.y < 0.0){
      float t = clamp(dot(p - inputVertices[2], inputVertices[0] - inputVertices[2]) / dot(inputVertices[0] - inputVertices[2], inputVertices[0] - inputVertices[2]), 0.0, 1.0);
      result[barycentricIndex] = vec3(t, 0.0, 1.0 - t);
    }else if(uvw.z < 0.0){
      float t = clamp(dot(p - inputVertices[0], inputVertices[1] - inputVertices[0]) / dot(inputVertices[1] - inputVertices[0], inputVertices[1] - inputVertices[0]), 0.0, 1.0);
      result[barycentricIndex] = vec3(1.0 - t, t, 0.0);
    }
  }

  // Calculate the output vertices by multiplying the barycentric coordinates with the input vertices.
  outputVertices[0] = (inputVertices[0] * result[0].x) + (inputVertices[1] * result[0].y) + (inputVertices[2] * result[0].z);
  outputVertices[1] = (inputVertices[0] * result[1].x) + (inputVertices[1] * result[1].y) + (inputVertices[2] * result[1].z);
  outputVertices[2] = (inputVertices[0] * result[2].x) + (inputVertices[1] * result[2].y) + (inputVertices[2] * result[2].z);

  // Return the barycentric coordinates for other possible calculations like texture coordinate interpolation and the like.
  return result;
}

// ...

  uint vertexIndex = uint(gl_VertexIndex);
  
  if(vertexIndex < countTotalVertices){ 

    const float GoldenRatio = 1.61803398874989485, // (1.0 + sqrt(5.0)) / 2.0 (golden ratio)
                IcosahedronLength = 1.902113032590307, // sqrt(sqr(1) + sqr(GoldenRatio))
                IcosahedronNorm = 0.5257311121191336, // 1.0 / IcosahedronLength
                IcosahedronNormGoldenRatio = 0.85065080835204; // GoldenRatio / IcosahedronLength

    const vec3 faceVertices[12] = vec3[12](
      vec3(0.0, IcosahedronNorm, IcosahedronNormGoldenRatio),
      vec3(0.0, -IcosahedronNorm, IcosahedronNormGoldenRatio),
      vec3(IcosahedronNorm, IcosahedronNormGoldenRatio, 0.0),
      vec3(-IcosahedronNorm, IcosahedronNormGoldenRatio, 0.0),
      vec3(IcosahedronNormGoldenRatio, 0.0, IcosahedronNorm),
      vec3(-IcosahedronNormGoldenRatio, 0.0, IcosahedronNorm),
      vec3(0.0, -IcosahedronNorm, -IcosahedronNormGoldenRatio),
      vec3(0.0, IcosahedronNorm, -IcosahedronNormGoldenRatio),
      vec3(-IcosahedronNorm, -IcosahedronNormGoldenRatio, 0.0),
      vec3(IcosahedronNorm, -IcosahedronNormGoldenRatio, 0.0),
      vec3(-IcosahedronNormGoldenRatio, 0.0, -IcosahedronNorm),
      vec3(IcosahedronNormGoldenRatio, 0.0, -IcosahedronNorm)
    );

    const uvec3 faceIndices[20] = uvec3[20](
      uvec3(0u, 5u, 1u), uvec3(0u, 3u, 5u), uvec3(0u, 2u, 3u), uvec3(0u, 4u, 2u), uvec3(0u, 1u, 4u),
      uvec3(1u, 5u, 8u), uvec3(5u, 3u, 10u), uvec3(3u, 2u, 7u), uvec3(2u, 4u, 11u), uvec3(4u, 1u, 9u),
      uvec3(7u, 11u, 6u), uvec3(11u, 9u, 6u), uvec3(9u, 8u, 6u), uvec3(8u, 10u, 6u), uvec3(10u, 7u, 6u),
      uvec3(2u, 11u, 7u), uvec3(4u, 9u, 11u), uvec3(1u, 8u, 9u), uvec3(5u, 10u, 8u), uvec3(3u, 7u, 10u)
    );

    uint triangleIndex = vertexIndex / 3u,   
         triangleVertexIndex = vertexIndex - (triangleIndex * 3u),         
         faceIndex = triangleIndex / squaredResolution,
         triangleSubdivisionIndex = triangleIndex - (faceIndex * squaredResolution);
         
    uvec3 faceVertexIndices = faceIndices[faceIndex];
    
    vec3 inputVertices[3] = vec3[3]( 
      faceVertices[faceVertexIndices.x], 
      faceVertices[faceVertexIndices.y], 
      faceVertices[faceVertexIndices.z]
    );

    vec3 outputVertices[3];

    tessellateTriangle(triangleSubdivisionIndex, resolution, inputVertices, outputVertices);

#ifdef CCW
    sphereNormal = normalize(outputVertices[2u - triangleVertexIndex]);
#else
    sphereNormal = normalize(outputVertices[triangleVertexIndex]);
#endif

    // 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));

// ...
    
   }

// ...