Mathematically producing sphere-shaped hexagonal grid

Question

I am trying to create a shape similar to this, hexagons with 12 pentagons, at an arbitrary size.

(Image Source)

The only thing is, I have absolutely no idea what kind of code would be needed to generate it!

The goal is to be able to take a point in 3D space and convert it to a position coordinate on the grid, or vice versa and take a grid position and get the relevant vertices for drawing the mesh.

I don't even know how one would store the grid positions for this. Does each "triagle section" between 3 pentagons get their own set of 2D coordinates?

I will most likely be using C# for this, but I am more interested in which algorithms to use for this and an explanation of how they would work, rather than someone just giving me a piece of code.

Answer 1

First some analysis of the image in the question: the spherical triangle spanned by neighbouring pentagon centers seems to be equilateral. When five equilateral triangles meet in one corner and cover the whole sphere, this can only be the configuration induced by a icosahedron. So there are 12 pentagons and 20 patches of a triangular cutout of a hexongal mesh mapped to the sphere.

So this is a way to construct such a hexagonal grid on the sphere:

1. Create triangular cutout of hexagonal grid: a fixed triangle (I chose (-0.5,0),(0.5,0),(0,sqrt(3)/2) ) gets superimposed a hexagonal grid with desired resolution n s.t. the triangle corners coincide with hexagon centers, see the examples for n = 0,1,2,20:

2. Compute corners of icosahedron and define the 20 triangular faces of it (see code below). The corners of the icosahedron define the centers of the pentagons, the faces of the icosahedron define the patches of the mapped hexagonal grids. (The icosahedron gives the finest regular division of the sphere surface into triangles, i.e. a division into congruent equilateral triangles. Other such divisions can be derived from a tetrahedron or an octahedron; then at the corners of the triangles one will have triangles or squares, resp. Furthermore the fewer and bigger triangles would make the inevitable distortion in any mapping of a planar mesh onto a curved surface more visible. So choosing the icosahedron as a basis for the triangular patches helps minimizing the distortion of the hexagons.)

3. Map triangular cutout of hexagonal grid to spherical triangles corresponding to icosaeder faces: a double-slerp based on barycentric coordinates does the trick. Below is an illustration of the mapping of a triangular cutout of a hexagonal grid with resolution n = 10 onto one spherical triangle (defined by one face of an icosaeder), and an illustration of mapping the grid onto all these spherical triangles covering the whole sphere (different colors for different mappings):

Here is Python code to generate the corners (coordinates) and triangles (point indices) of an icosahedron:

from math import sin,cos,acos,sqrt,pi
s,c = 2/sqrt(5),1/sqrt(5)
topPoints = [(0,0,1)] + [(s*cos(i*2*pi/5.), s*sin(i*2*pi/5.), c) for i in range(5)]
bottomPoints = [(-x,y,-z) for (x,y,z) in topPoints]
icoPoints = topPoints + bottomPoints
icoTriangs = [(0,i+1,(i+1)%5+1) for i in range(5)] +\
             [(6,i+7,(i+1)%5+7) for i in range(5)] +\
             [(i+1,(i+1)%5+1,(7-i)%5+7) for i in range(5)] +\
             [(i+1,(7-i)%5+7,(8-i)%5+7) for i in range(5)]

And here is the Python code to map (points of) the fixed triangle to a spherical triangle using a double slerp:

# barycentric coords for triangle (-0.5,0),(0.5,0),(0,sqrt(3)/2)
def barycentricCoords(p):
  x,y = p
  # l3*sqrt(3)/2 = y
  l3 = y*2./sqrt(3.)
  # l1 + l2 + l3 = 1
  # 0.5*(l2 - l1) = x
  l2 = x + 0.5*(1 - l3)
  l1 = 1 - l2 - l3
  return l1,l2,l3

from math import atan2
def scalProd(p1,p2):
  return sum([p1[i]*p2[i] for i in range(len(p1))])
# uniform interpolation of arc defined by p0, p1 (around origin)
# t=0 -> p0, t=1 -> p1
def slerp(p0,p1,t):
  assert abs(scalProd(p0,p0) - scalProd(p1,p1)) < 1e-7
  ang0Cos = scalProd(p0,p1)/scalProd(p0,p0)
  ang0Sin = sqrt(1 - ang0Cos*ang0Cos)
  ang0 = atan2(ang0Sin,ang0Cos)
  l0 = sin((1-t)*ang0)
  l1 = sin(t    *ang0)
  return tuple([(l0*p0[i] + l1*p1[i])/ang0Sin for i in range(len(p0))])

# map 2D point p to spherical triangle s1,s2,s3 (3D vectors of equal length)
def mapGridpoint2Sphere(p,s1,s2,s3):
  l1,l2,l3 = barycentricCoords(p)
  if abs(l3-1) < 1e-10: return s3
  l2s = l2/(l1+l2)
  p12 = slerp(s1,s2,l2s)
  return slerp(p12,s3,l3)