snowflake_fractal.py

# Imports
import matplotlib.pyplot as plt
from matplotlib.path import Path
import matplotlib.patches as patches
import numpy as np
import math
import time

# Define Constants
pi = math.pi

codes = [
    Path.MOVETO,
    Path.LINETO,
    Path.LINETO,
	Path.CLOSEPOLY
]



# Methods for Constructing Triangles

"""
* Function build_rotation_matrix
* Input:	theta <float> - the radians by which the vector should rotate
* Out:		<np.array(2,2)> - the rotation matrix given theta
"""
def build_rotation_matrix(theta):
	s, c = np.sin(theta), np.cos(theta)
	return np.array([[c, -s],
					 [s,  c]])

"""
* Function build_radius_vector
* Input:	a <float> - the length of the radius
* Out		<np.array(2,1)> - the radius vector
"""
def build_radius_vector(a):
	x = a * math.sqrt(3) / 3
	return np.array([x, 0])

"""
* Function rotate_radius
* Input:	radius <np.array> - the radius vector
* 			rotation <np.array> - the rotation matrix
* Out:		<np.array> - the rotated radius vector
"""
def rotate_radius(radius, rotation):
	return np.dot(rotation, radius)

"""
* Function build_triangle
* Input:	center <np.array> - specifies the center of triangle
* 			a <float> - the length of the sides of the triangle
*			theta_null <float> - this is the offset from the x axis of
								 the angle the "center" vector makes
* Out:		List<np.array> - a list of points for the triangle
"""
def build_triangle(center, a, theta_null=0):

	# Define Radius Rotations
	thetas = [
		pi / 2,
		(7 * pi) / 6,
		(11 * pi) / 6
	]
	# Build Radius Phasor
	radius = build_radius_vector(a)

	# Rotate the Radius Phasor to build the points
	rotations = [build_rotation_matrix(theta_null + theta) for theta in thetas]
	radian = [rotate_radius(radius, rotation) for rotation in rotations]

	return [center + radius for radius in radian] + [np.array([0,0])]

"""
* Function get_next_centers
* Input:	center <np.array> - specifies the center of triangle
* 			a <float> - the length of the sides of the triangle
*			theta_null <float> - this is the offset from the x axis of
								 the angle the "center" vector makes
* Out:		List<np.array> - a list of points for the next 3 centers
"""
def get_next_centers(center, a, theta_null=0):

	r1 = (a * math.sqrt(3)) / 6
	r2 = r1 + r1 /2
	r = r1 + r2

	thetas = [
		pi / 6,
		(5 * pi) / 6,
		(3 * pi) / 2
	]

	radius = build_radius_vector(r)
	rotations = [build_rotation_matrix(theta_null + theta) for theta in thetas]
	radian = [rotate_radius(radius, rotation) for rotation in rotations]
	return [center + radius for radius in radian]

"""
* Function assemble_path
* Input:	points List<np.array> - a list of points from which to build a path
* Out:		<Path> - a matplotlib.path object
"""
def assemble_path(points):
	return Path(points, codes)

"""
* Function graph
* Input:	List<Path> - list of path objects
* Out:		None
"""
def graph(paths):
	fig = plt.figure()
	ax = fig.add_subplot(111)

	for path in paths:
		patch = patches.PathPatch(path, facecolor='black', lw=0, ec=None)
		ax.add_patch(patch)

	ax.set_xlim(-2,2)
	ax.set_ylim(-2,2)
	plt.show()

if __name__ == '__main__':
	
	centers = [np.array([0,0])]
	a = 2.0
	paths = []

	k = 0
	start = time.time()
	while k < 12:
		iter_time = time.time()
		new_centers = []
		for center in centers:
			theta_null = 0 if k % 2 == 0 else (-pi / 3)
			points = build_triangle(center, a, theta_null)
			paths.append(assemble_path(points))
			new_centers += get_next_centers(center, a, theta_null)
		print("Layer %d [%0.5fs]" % (k, time.time() - iter_time))
		centers = new_centers
		a /= 2
		k += 1
	print("Finish [%0.2fs]" % (time.time() - start))
		
	graph(paths)