# New addon: Snowflake generator

A post by a Google+ contact of mine recently caught my eye with a cool animated gif.

It showed a snowflake fractal rendered in 2D with growing complexity. I love fractals (who doesn’t?), and snowflake shapes are among the most beautiful of them. The post said that the image was generated by a nice little web-application, developed by Italian computer science student, Marco Cameriero.

His app generates what’s called a “Koch curve” or “Koch Snowflake” shape. I was intrigued and read more about these shapes in Wikipedia. The procedure to generate such a shape seemed fairly straightforward, so I decided to try and write a blender script that generates Koch Snowflake meshes.

The problem is that I’m not exactly an expert in trigonometry or linear algebra, so calculating the equilateral triangle based vertex positions was a bit of a handful for me.

Luckily, blender’s already has objects that take care of all the math behind the scenes (vectors, matrices, etc), so I found ways to use blender’s transformation tools and some very basic vector operations to avoid calculating vertex positions myself.

I’ll explain more about the algorithm in a moment, but first, for those not interested in Blender Python programming, here’s a download link for the Addon: DOWNLOAD.

For detailed installation and usage, see the WIKI PAGE.

**How to use the script:**

**Updated: Video tutorial now added**

**Building the snowflake in Python:**

To create a Koch snowflake, you need:

1. A basic polygonal mesh (triangle, rectangle, pentagon, etc ).

2. An object that will link this mesh geometry to the scene.

3. To iterate over the polygon’s edges, split them and create triangles at the center of each edge.

4. Repeat this process for as many iterations as specified by the user.

Steps 1 and 2 were solved before and all I had to do is reuse code written by others (open source rulez!!):

The calculations necessary to create a basic polygon already appear in the “curvaceous galore” add-on that comes bundled with Blender (thank you Jimmy Hazevoet & testscreenings!). The function I borrowed code from to do this is “StarCurve” (line 258).

The code for creating a new mesh from a provided list of vertices and edges, and creating an object that links this geometry to the scene is specified in the awesome Blender 2.5 API “Code Snippets” (thanks to the authors!).

**Algorithm for creating the Koch snowflake from a basic polygon:**

Since I’m pretty lame in math, I shied away from trying to calculate the entire shape, and the positions of each vertex based on equilateral triangle trigonometry. Instead, I went for an automation algorithm for creating this shape, that would mimic (to some degree at least) a manual modeling process.

The steps required to create the snowflake from a polygon are these:

1. Iterate over each edge.

2. Split each edge to 3 equal edges (using the mesh.subdivide operator).

3. Select the edge at the center of the 3 newly created edges.

4. Subdivide it again to two equal edges.

5. Select the central vertex.

6. Translate this vertex to a position and distance that will create an equilateral triangle between it and the other two vertices of the subdivided edge. This does require some calculations, which I’ll discuss later.

These 6 steps (if performed on all edges) define one “iteration” of the whole snowflake generating algorithm. The next iteration will go over all the new edges and do the same with each, thus creating a more complex fractal shape.

**Concerns and possible caveats:**

I had one major concern regarding this approach: when you subdivide an edge, you create new edges. How does this affect the “edges” collection? Are the new edges simply added at the end of the line, or is the whole list re-arranged?

And if we create a reference object to all mesh edges ( for instance: “edges = bpy.context.object.data.edges” ), will that reference still be valid after the subdivision, or will it become stale?

And what about the vertices? An edge is simply an index that’s associated with the pair of vertices (or vertex indices) that define it. So what if the vertices’ indices are rearranged when new verts are added?

It was pretty clear I can’t simply iterate over the mesh’s edges (object.data.edges), since this list will change dynamically with subdivisions.

To solve this, I created a static dictionary, that has an arbitrary edge index as a key, and a list containing the associated vertex indices as the value.

Example: { 0 : [ 0, 1 ], 1 : [ 1, 2 ], 2: [2, 0 ] }

But the concerns mentioned above are just as relevant to vertices as they are to edges: what happens when a new vertex is added due to a subdivision? is it added at the end of the list, or are all indices recalculated?

I tested it quickly in the Python Console, and found out that new vertices are added at the end, regardless of their position on the mesh. Yay!

If this wasn’t the case, I’d have no way to store the undivided edges’ information and iterate only over them in each fractal iteration. Which means I’d have no choice but to create all the verts manually and calculate their positions with complex math (unlikely).

**Finding central edges and vertices via Set boolean operations**:

After you subdivide an edge into 3 edges, all three are selected. You have no immediate simple way to know which one is the central edge, which must be subdivided again.

One way I came up with to identify the central edge, is to compare the vertex indices between the three edges. The central edge shares both its verts with the other two edges, while the peripheral edges share only one vertex with another edge. So if I find the edge with two shared vertices, I’ll find the central edge.

A relatively quick way to tell how many shared verts an edge has, is by using Set boolean operations.

If you have two sets, you can easily perform an intersection, union or symmetrical difference between the sets to get the shared elements, all elements or only the non-shared elements, accordingly:

o = bpy.context.object # Selected object d = o.data # Object's data (geometry) bm = bmesh.from_edit_mesh( d ) # Create BMESH object from object data bv = bm.verts # Reference object's verts specifically ## Find innermost edge # Create dictionary of selected edges' indices and their verts selected_edges = { e.index : set([ v.index for v in e.verts ]) \ for e in bm.edges if e.select } # Innermost edge shares both its verts with the other edges mid = 42 # Initialize with out-of-range index to to elicit error of algorithm fails for i in selected_edges.keys(): other_edges = set( selected_edges.keys() ) - set( [i] ) # Set Symmetrical Difference # check how many joint verts this edge has with the others num_of_joint_verts = 0 for oe in other_edges: if selected_edges[i] & selected_edges[oe]: # Intersect vertex sets num_of_joint_verts += 1 # The middle edge will have two joint verts if num_of_joint_verts == 2: mid = i # Record middle edge's index break # Select innermost edge (by selecting its verts) bpy.ops.mesh.select_all(action = 'DESELECT') # Deselect all verts for v in selected_edges[ mid ]: bv[ v ].select = True bm.select_flush( True )

I used a similar technique to find the innermost vertex created after the central edge is subdivided again.

The edge is split into two edges, that have 3 vertices. The two edges share only one vertex, which is the central vertex.

So an intersection between the vertex sets of both edges will discover the central vertex:

## Select innermost vert (shared vert between both edges) # Create list of sets, each containing the indices of a selected edge's verts selected_edges_verts = [ set([ v.index for v in e.verts ]) for e in bm.edges if e.select ] # Find common vert by intersecting both vert sets then poping out the only member of the set joint_vert = list( selected_edges_verts[0] & selected_edges_verts[1] ).pop() # Find the verts that aren't joint (via set symmetrical_difference) other_verts = list( selected_edges_verts[0] ^ selected_edges_verts[1] )

**Translating the central vertex to create an equilateral triangle:
** At first I did try to figure out the trigonometric relationships between the edges of an equilateral triangle, and maybe find a way to calculate the new position of the central vertex. But I didn’t find a quick solution (probably due to my lack of knowledge), and decided to try a different solution to this problem.

What I needed was in essence a translation vector. This needs to point outwards from the object origin and carry the vertex the distance require to create an equilateral triangle. To get that vector, first I subtracted the location vectors (coordinates) of the two peripheral vertices.

This resulted in a vector with the required length, that lies perpendicular to the Vector I needed to translate the vertex.

I then multiplied it with a rotation matrix (that can be created rather simply using a Matrix object), to change the direction of the Vector.

This resulted in a translation vector that I could use with the transform.translate operator:

# Calculate its new position: should create an equilateral triangle vdiff = bv[ other_verts[1] ].co - bv[ other_verts[0] ].co vrot = vdiff * rot_matrix # Select middle vert and translate it by the rotate vector bpy.ops.mesh.select_all(action = 'DESELECT') bv[ joint_vert ].select = True bm.select_flush( True ) bpy.ops.transform.translate( value = tuple( vrot ) )

There’s another small piece of code I didn’t show here, meant to make sure the triangle’s edge is directed outwards and not inwards, but I think I probably digressed enough about the code already, and you’re welcome to explore the source further to see the other tricks I used.

**Creating more complex models like the one displayed at the top of this post:**

Obviously, there’s only so much you can do with an edge loop that looks like a snowflake fractal’s outline.

I might expand on this script one day to create faceted shapes and profiles, but at the moment I’m pretty satisfied with it, as it provides a basis to work with.

The shape shown at the featured image above has obviously been processed beyond what the script provides. Here’s how I did it:

1. After creating a hexagon based fractal with two iterations, I selected the all the verts and duplicated them.

2. I then scaled these duplicated verts inwards, then used the “Smooth Vertex” operator in them (can be found in the Tools Panel under “Deform”, or in the Special’s menu (shortcut: “W”) –> Smooth.

3. After getting a nice inner hub I used the Skin Modifier to create a 3D shape out of these two shapes.

4. I reduced the radius of each vertex using the “Ctrl + A” shortcut with all verts selected, and then applied the skin modifier.

5. The only thing left was to Bridge (loop tools) between areas of the inner and outer hubs (I chose two faces from each hub and bridged between them). To create a tapered bridge I set the bridge segments to 3 and the Strength to around 0.5 in Cubic interpolation mode.

Another probably faster way to move that central vertice will be this:

After subtracting the two peripheral vertices and getting the required length, L. You find the normal of the central vertices. The Blender Python API should have a function for that (I am not good with the API).

Then u multiply the lenght, L, you got above with sqrt(3)/2 to get the magnitude of lenght u should move the central vertex along its normal. Lets call this M.

U then multipy M with a unit vector along the normal of that central vertex. And move the central vertex along the vector u get.

This is assuming the Blender Python API has a function for finding the normals of vertices. I know Blender can show the normal of vertices in the 3D viewport.

Nice add-on by the way. I am fascinated by fractals. Thanks for the add-on.

Cheers 🙂

Thanks for the tip, Phillip – I’ll check it out.

At some point I’ll probably want to add more features to this script (such as adding more complex polygons to edges, and some other cool things Marco Cameriero did in his web-app), and this might prove to be a better method.

Cheers!