Meshes with Python & Blender: Icospheres

In the third part of the series we get into making icosahedra, subdividing and refining them to spheres. To top it off, we’ll also look at two ways of setting a mesh’s shading to smooth.

Tutorial Series

What is a icosahedron

An icosahedron is a polyhedron with 20 faces. There are several kinds of icosahedra. However, to make icospheres we’ll only be looking at convex regular icosahedra (also the most famous kind). You can find more about them and their properties on Wikipedia

So, why icospheres? Icospheres have a more even distribution of geometry than UV Spheres. Deforming UV Spheres often gives strange results near the poles due to the higher density of geometry, while icospheres give a more even and organic looking result. On top of that icospheres are asymmetrical which helps sell an organic deformation.

This tutorial is based on the original icosahedron code from Andreas Kahler adapted to Python 3 and Blender.

Setting it up

I bet you already know this by now. Let’s start importing and then move on to our usual scaffolding.

import bpy
from math import sqrt

# -----------------------------------------------------------------------------
# Settings
name = 'Icosomething'
scale = 1
subdiv = 5

# -----------------------------------------------------------------------------
# Add Object to Scene

mesh = bpy.data.meshes.new(name)
mesh.from_pydata(verts, [], faces)

obj = bpy.data.objects.new(name, mesh)
bpy.context.scene.objects.link(obj)

bpy.context.scene.objects.active = obj
obj.select = True

In the settings, subdiv will control how many times to subdivide the mesh and scale will be a simple uniform scale parameter much like the ones in the previous tutorial. Setting subdiv to zero will create a icosahedron (instead of an icosphere). Note that a subdiv value of 9 will result in a mesh with over 5 million faces. You may want to stay below that value depending on your hardware (or how badly you want to crash Blender!).

Putting the sphere in icosphere

Simply subdividing the icosahedron will only get us a refined icosahedron. We need to make sure the vertices come together in a way that resembles a sphere.

To make this happen we need to make sure the vertices we add lie on the unit sphere. The unit sphere is an “imaginary” sphere of radius 1. We can determine the position of each point (vertex) on the unit sphere with a simple formula and then fix the coordinates to it. More on Wikipedia

To do this we’ll have a vertex() function that fix to the unit sphere (and does the scaling too).

def vertex(x, y, z):
    """ Return vertex coordinates fixed to the unit sphere """

    length = sqrt(x**2 + y**2 + z**2)

    return [(i * scale) / length for i in (x,y,z)]

Make the base Icosahedron

Now that we know the verts are falling on the unit sphere we can go ahead and create the base icosahedron. Like the cube before, the easiest way is to input the vertices and faces manually.

One of the ways to build an icosahedron is to consider its vertices as the corners of three orthogonal golden planes. They are golden planes because they follow the golden ratio. The vertices of these planes lie on the coordinates (0, ±1, ±φ), (±φ, 0, ±1) and (±1, ±φ, 0). Note that the letter φ (phi) represents the golden ratio value, while ± means “negative or positive”.

These combinations result in 12 vertices, which create 20 equilateral triangles with 5 triangles meeting at each vertex. Check out the following diagram (points are numerated for the faces list below).

Damn, that’s a lot of math! However it’s all fairly straight forward (and repetitive) when you think about it. A math jock would probably find this boring. Here’s the code to make the icosahedron.

# --------------------------------------------------------------
# Make the base icosahedron

# Golden ratio
PHI = (1 + sqrt(5)) / 2

verts = [
          vertex(-1,  PHI, 0),
          vertex( 1,  PHI, 0),
          vertex(-1, -PHI, 0),
          vertex( 1, -PHI, 0),

          vertex(0, -1, PHI),
          vertex(0,  1, PHI),
          vertex(0, -1, -PHI),
          vertex(0,  1, -PHI),

          vertex( PHI, 0, -1),
          vertex( PHI, 0,  1),
          vertex(-PHI, 0, -1),
          vertex(-PHI, 0,  1),
        ]


faces = [
         # 5 faces around point 0
         [0, 11, 5],
         [0, 5, 1],
         [0, 1, 7],
         [0, 7, 10],
         [0, 10, 11],

         # Adjacent faces
         [1, 5, 9],
         [5, 11, 4],
         [11, 10, 2],
         [10, 7, 6],
         [7, 1, 8],

         # 5 faces around 3
         [3, 9, 4],
         [3, 4, 2],
         [3, 2, 6],
         [3, 6, 8],
         [3, 8, 9],

         # Adjacent faces
         [4, 9, 5],
         [2, 4, 11],
         [6, 2, 10],
         [8, 6, 7],
         [9, 8, 1],
]

Strategy for subdivision

We can grab a triangle and split each edge creating three triangles in its place. Basically turning triangles into little triforces. Note that when I say split I’m not talking about actually running an operator and cutting the edge, but rather placing a new vertex in the middle of the other original two and making three new faces.

However if we went around splitting all edges we would quickly run into the same edges we have already split. This would result in a lot of duplicated verts and a headache when trying to build the faces.

To prevent that, let’s keep a list of the edges we have already split (a cache), and check it before splitting. This cache will be a dictionary. The keys will be the index of the vertices, ordered from smaller to greater. That way the key will remain the same, no matter how we pass the edge’s vertices.

middle_point_cache = {}

def middle_point(point_1, point_2):
    """ Find a middle point and project to the unit sphere """

    # We check if we have already cut this edge first
    # to avoid duplicated verts
    smaller_index = min(point_1, point_2)
    greater_index = max(point_1, point_2)

    key = '{0}-{1}'.format(smaller_index, greater_index)

    if key in middle_point_cache:
        return middle_point_cache[key]

    # If it's not in cache, then we can cut it
    vert_1 = verts[point_1]
    vert_2 = verts[point_2]
    middle = [sum(i)/2 for i in zip(vert_1, vert_2)]

    verts.append(vertex(*middle))

    index = len(verts) - 1
    middle_point_cache[key] = index

    return index

The middle vertex is calculated by adding the coordinates from both vertices and dividing them by 2. Finally we add it to the cache and return the index to make the faces list.

Subdividing

With the middle_point() function done we can move over to looping and making the subdivisions.

At each subdivisions level we create a new empty list of faces, and at the end we replace the original faces list with the new one. Then we go through each face, finding the middle point for it’s three edges, storing the indices and creating 4 new faces from them (remember the diagram above).

# Subdivisions
# --------------------------------------------------------------

for i in range(subdiv):
    faces_subdiv = []

    for tri in faces:
        v1 = middle_point(tri[0], tri[1])
        v2 = middle_point(tri[1], tri[2])
        v3 = middle_point(tri[2], tri[0])

        faces_subdiv.append([tri[0], v1, v3])
        faces_subdiv.append([tri[1], v2, v1])
        faces_subdiv.append([tri[2], v3, v2])
        faces_subdiv.append([v1, v2, v3])

    faces = faces_subdiv

Make it smooth

We now have a mesh that approximates a sphere quite well, but still looks faceted. Time to smooth it out.

Smooth shading is an attribute of faces. So, to set a whole mesh as smooth we need to enable smooth shading on all of it’s faces. You can do this using the same operator we use when we click the Set smooth button in the UI:

bpy.ops.object.shade_smooth()

This will work just fine for this script because the context is right for the operator. In other cases, you might find the operator refuses to run due to an “incorrect context”. context in Blender is a sort of god-variable that contains information about the current state of the application. This includes things like mouse cursor position, the active area and more. You can override the context when calling an operator but currently there’s no easy way of knowing what each operator wants to see in the context.

Luckily there’s another way to do this the “low-level” way by setting each face in a loop.

for face in mesh.polygons:
    face.use_smooth = True

In the world of Blender scripting “Low Level” refers to skipping operators and going straight to the methods and attributes of objects. from_pydata() is another example of low-level work.

The benefits of low-level is that it doesn’t depend on context, it’s often more flexible and saves you the overhead of going through the operators system. In this case you could also decide to apply smooth only on some faces.

Final Code

import bpy
from math import sqrt

# -----------------------------------------------------------------------------
# Settings

scale = 1
subdiv = 5
name = 'Icosomething'


# -----------------------------------------------------------------------------
# Functions

middle_point_cache = {}


def vertex(x, y, z):
    """ Return vertex coordinates fixed to the unit sphere """

    length = sqrt(x**2 + y**2 + z**2)

    return [(i * scale) / length for i in (x,y,z)]


def middle_point(point_1, point_2):
    """ Find a middle point and project to the unit sphere """

    # We check if we have already cut this edge first
    # to avoid duplicated verts
    smaller_index = min(point_1, point_2)
    greater_index = max(point_1, point_2)

    key = '{0}-{1}'.format(smaller_index, greater_index)

    if key in middle_point_cache:
        return middle_point_cache[key]

    # If it's not in cache, then we can cut it
    vert_1 = verts[point_1]
    vert_2 = verts[point_2]
    middle = [sum(i)/2 for i in zip(vert_1, vert_2)]

    verts.append(vertex(*middle))

    index = len(verts) - 1
    middle_point_cache[key] = index

    return index


# -----------------------------------------------------------------------------
# Make the base icosahedron

# Golden ratio
PHI = (1 + sqrt(5)) / 2

verts = [
          vertex(-1,  PHI, 0),
          vertex( 1,  PHI, 0),
          vertex(-1, -PHI, 0),
          vertex( 1, -PHI, 0),

          vertex(0, -1, PHI),
          vertex(0,  1, PHI),
          vertex(0, -1, -PHI),
          vertex(0,  1, -PHI),

          vertex( PHI, 0, -1),
          vertex( PHI, 0,  1),
          vertex(-PHI, 0, -1),
          vertex(-PHI, 0,  1),
        ]


faces = [
         # 5 faces around point 0
         [0, 11, 5],
         [0, 5, 1],
         [0, 1, 7],
         [0, 7, 10],
         [0, 10, 11],

         # Adjacent faces
         [1, 5, 9],
         [5, 11, 4],
         [11, 10, 2],
         [10, 7, 6],
         [7, 1, 8],

         # 5 faces around 3
         [3, 9, 4],
         [3, 4, 2],
         [3, 2, 6],
         [3, 6, 8],
         [3, 8, 9],

         # Adjacent faces
         [4, 9, 5],
         [2, 4, 11],
         [6, 2, 10],
         [8, 6, 7],
         [9, 8, 1],
        ]


# -----------------------------------------------------------------------------
# Subdivisions

for i in range(subdiv):
    faces_subdiv = []

    for tri in faces:
        v1 = middle_point(tri[0], tri[1])
        v2 = middle_point(tri[1], tri[2])
        v3 = middle_point(tri[2], tri[0])

        faces_subdiv.append([tri[0], v1, v3])
        faces_subdiv.append([tri[1], v2, v1])
        faces_subdiv.append([tri[2], v3, v2])
        faces_subdiv.append([v1, v2, v3])

    faces = faces_subdiv


# -----------------------------------------------------------------------------
# Add Object to Scene

mesh = bpy.data.meshes.new(name)
mesh.from_pydata(verts, [], faces)

obj = bpy.data.objects.new(name, mesh)
bpy.context.scene.objects.link(obj)

bpy.context.scene.objects.active = obj
obj.select = True


# -----------------------------------------------------------------------------
# Smoothing

#bpy.ops.object.shade_smooth()

for face in mesh.polygons:
    face.use_smooth = True

Wrap up

That completes the third part of this series. If you are interested in making spheres or apply points/objects in a spherical way consider reading more about the unit sphere and how to apply it along normals.

Things you can do for yourself:

  • Optimize it (Hint: you don’t need to store the key as a string)
  • Bring the rotation and translation matrices from previous part
  • Remove the scale setting and use a matrix instead

Next time we’ll go back to cubes and learn how to make a round cube, apply modifiers and more. Do you have any questions, or suggestions for the next tutorials? Leave a comment below!