Drawing the Koch Snowflake

The Koch Snowflake is a fractal defined by repeatedly adding spikes to the edges of an equilateral triangle. The middle third of each edge is doubled and bent outward, forming a smaller equilateral triangle pointing out of the edge. The length of the curve grows by a factor of 4/3 longer.

Repeat that on the new edges, and then on their new edges, forever. Each time, the curve grows longer by a factor 4/3, resulting in an infinite final length! The Hausdorff dimension of the Koch curve is log₃(4) = 1.2618595...

An Iterated Equation

The Koch Snowflake is typically drawn recursively. An array of edges is subdivided, exactly as described in the definition of the fractal, starting with three edges of an equilateral triangle. However, this lacks information about the interior and exterior of the curve. The Mandelbrot set and other fractals are drawn iteratively, usually by coloring points by how fast they diverge toward infinity (escape time). This reveals beautiful patterns and structure in the area near the boundary, just outside. Here is a graph of the Koch Snowflake:

My approach is to construct an iterated equation for a set whose boundary is the Koch Snowflake curve, using hexagonal symmetry. Because the fractal is constructed from equilateral triangles, and two equilateral triangles overlaid form a hexagon, each branch of the flake also exhibits hexagonal symmetry. This converges twice as fast as using triangles.

To iterate a point, Z, it tests if it's inside the middle hexagon. If it is, then it's inside the set, and it's done. Otherwise, it checks to see which of the six smaller hexagons it's closest to, and then reflects and scales the point to place that smaller hexagon at the same size and position as the largest one. This process is iterated until either the point is finally transformed into the largest hexagon, meaning it's inside the fractal set, or it runs for the maximum specified number of iterations, and is probably outside the set.

To simplify this approach (at a small performance cost), the point can be repeatedly rotated by an angle of ± π3, and then only tested against the small hexagon on the right when it comes around into position. Stated precisely, the reflection must happen when a point's angle is |θ| < π6, and it's x coordinate is greater than 1. This can be written with complex numbers simply as Re(z³) > 0 and Re(z) > 1.

In polar notation, z³ = (e^iθ)³ = e^3iθ

Therefore, θ = ±π6 ⇒ 3θ = ±π2, which is just the vertical axis, so testing if Re(z³) > 0 and Re(z) > 0 is equivalent to testing if |θ| < π6, but much faster to compute.

Finally, since the reflection must happen only if the point is outside the center hexagon, whose right edge is at 1, it must also check if Re(z) > 1. This simplifies further when combined with the previous condition, since Re(z) > 0 disappears, being weaker than Re(z) > 1

Therefore, the condition for reflection is Re(z³) > 0 and Re(z) > 1, which is true in the red region shown below. It covers only the small hexagon on the right, and its descendents, which covers any possible path into smaller hexagons in future iterations.

Finally, the reflection operation itself is just scaling and translation. It must invert the fractal operation to map the smaller hexagon to the largest one. Recall in the definition of the fractal, the middle third of each edge becomes the new edge for the next iteration. Therefore, the inverse must scale by 3.

The center position of the small hexagon is 1 + 13, which after scaling by 3 is 3 + 33 = 4. Therefore, the transformation 4 − 3z maps it to the larger hexagon.

Putting it all together as an iterated formula results in the following code. This draws the interior of the snowflake set, where the value of z is evaluated for every pixel.

loop (n, 0, 64) {
if (re(z^3) > 0 && re(z) > 1) {
z = 4 - 3*z; // Reflect and scale.
}
z = z * e^(i*pi/3); // Rotate by π/3
}

Plotting the Escape Time

The iterated equation was derived for points inside the set, without considering the transformed value of z if the point is outside. However, since points outside the snowflake set are scaled in the wrong direction, they diverge toward infinity, which enables the formula to graph the escape time, exactly the same as for the Mandelbrot set.

The trick is to count iterations only when a reflection occurs, since the rotations do not change the magnitude of z.

complex k = 0; // Iteration counter.

loop (n, 0, 64) {
// Graph the escape time (k).
if (|z| > 60) return (k/9)^2*(z*0.004+1);

if (re(z^3) > 0 && re(z) > 1) {
z = 4 - 3*z; // Reflect and scale.
k = k+1;
}
z = z * e^(i*pi/3); // Rotate
}

return 0; // Draw black inside the set. ↗ Open the code in the Complex Plane Grapher.