Drawing the Koch snowflake!

.gitignore

*.pyc
.DS_Store

koch.py

#!/usr/bin/env python
# encoding: utf-8
from __future__ import division
import sys
import os
import Image, ImageDraw
import math
from collections import deque

"""
	Let's draw a fractal!
	
	How about the Koch snowflake? It seems simple and kind of neat.
	http://en.wikipedia.org/wiki/Koch_snowflake
"""

background = (0, 0, 0)

def segmentColor(age, maxAge,
                 youngColor = (255, 255, 255),
                 oldColor   = (128, 128, 128)):
	"""Returns the color to be used for a line of the given age."""
	
	if maxAge:
		ratio = age / maxAge
		return tuple(int((1 - ratio) * y + ratio * o) for (y, o) in zip(youngColor, oldColor))
	else:
		return youngColor
	

class Point(object):
	"""A point in two dimensions."""
	
	def __init__(self, x, y):
		self.x = x
		self.y = y
	
	def __repr__(self):
		return "%s(%.3f, %.3f)" % (type(self).__name__, self.x, self.y)
	
	def __add__(self, other):
		return Point(self.x + other.x, self.y + other.y)
	
	def __sub__(self, other):
		return Point(self.x / other.x, self.y / other.y)
	
	def __mul__(self, factor):
		return Point(self.x * factor, self.y * factor)
	
	def __truediv__(self, factor):
		return Point(self.x / factor, self.y / factor)
	
	def distanceFrom(self, other):
		return ((self.x - other.x) ** 2 + (self.y - other.y) ** 2) ** .5
	

class Segment(object):
	"""A line segment comprising two points."""
	
	def __init__(self, start, end, age=0):
		self.start = start
		self.end = end
		self.age = age
	
	def __str__(self):
		return "<%s of age %i at 0x%x>" % (type(self).__name__, self.age, id(self))
	

class KochSnowflake(object):
	"""A representation of an iteration of a simple fractal pattern."""
	
	def __init__(self, segments = None):
		if segments:
			self.segments = list(segments)
		else:
			side = 0.8 # large but not overflowing
			
			bottomLeft  = Point(-side, 0)
			bottomRight = Point(+side, 0)
			top = equalaterize(bottomRight, bottomLeft)
			
			offset = Point(0, -top.y / 3)
			# Let's center it in the square.
			
			bottomLeft  += offset
			bottomRight += offset
			top         += offset
			
			self.segments = [
				Segment(bottomRight, top),
				Segment(top, bottomLeft),
				Segment(bottomLeft, bottomRight)
			]
		
		self.age = max(segment.age for segment in self.segments)
	
	def advance(self, n = 1):
		"""Advance the fractal n iterations."""
		
		self.age += n
		
		for n in xrange(n):
			newSegments = deque()
			
			for segment in self.segments:
				points = [segment.end * (p / 3) + segment.start * ((3 - p) / 3) for p in xrange(4)]
				pointOut = equalaterize(points[1], points[2])
				
				newSegments.append(Segment(points[0], points[1], segment.age + 1))
				newSegments.append(Segment(points[1], pointOut))
				newSegments.append(Segment(pointOut, points[2]))
				newSegments.append(Segment(points[2], points[3], segment.age + 1))
			
			self.segments = newSegments
		
		return self # Chaining is fun!
	
	def render(self, size):
		"""
			Returns a picture of the fractal, bounded at ±1.0.
		"""
		
		image = Image.new("RGB", (size, size), background)
		drawLine = ImageDraw.Draw(image).line
		
		for segment in self.segments:
			drawLine(adaptPoint(segment.start, size) + adaptPoint(segment.end, size),
			         fill=segmentColor(segment.age, self.age))
		
		return image
	

def adaptPoint(point, size):
	"""Returns a 2-tuple of where a Point() bounded by ±1.0 would
	be located in an image with sides of size."""
	
	half = size // 2
	
	return tuple(half + half * v for v in (point.x, -point.y))

def equalaterize(a, b):
	"""Given two points return the third point
	needed to form an equalaterial triangle.
	
	Two possible points fit this criteria, returned is
	the one that would be above the two points if they
	were drawn on a vertical line with a to the left."""
	
	distance = a.distanceFrom(b)
	
	dX = b.x - a.x
	dY = b.y - a.y
	
	if dX:
		lineSlope = (b.y - a.y) / (b.x - a.x)
		lineAngle = math.atan(lineSlope)
	else:
		lineAngle = fullRotation / 2
	
	fullRotation = 2 * math.pi
	thirdRotation = fullRotation / 3
	
	# Vertical lines aren't handled. Expect problems.
	
	angleA = (lineAngle - thirdRotation) % fullRotation
	if angleA: slopeA = math.tan(angleA)
	
	angleB = (lineAngle + thirdRotation) % fullRotation
	if angleB: slopeB = math.tan(angleB)
	
	x = (b.y - a.y - b.x * slopeB + a.x * slopeA) / (slopeA - slopeB)
	y = a.y + (x - a.x) * slopeA
	
	return Point(x, y)
	

def main(*args):
	size = 800
	iterations = 5
	
	snowflake = KochSnowflake([Segment(Point(-1, 0), Point(+1, 0))])
	
	image = KochSnowflake().advance(iterations).render(size)
	
	image.save(os.path.expanduser("~/Desktop/test.png"), "PNG")

if __name__ == "__main__": sys.exit(main(*sys.argv[1:]))

mandelbrot.py

#!/usr/bin/env python
# encoding: utf-8
from __future__ import division, with_statement
import sys
import os
import Image
import time
from subprocess import Popen as Subprocess, PIPE

def colourer(intensity):
	"""Returns a nice colour given an intensity from 0 to 1."""
	
	try:
		return(colourer.cache[intensity])
		# Yes, this actually matters.
	except KeyError:
		r, g, b = 0, 0, 0
	
		if intensity:
			b = int(min(255 * (intensity * 3), 255))
		
			if intensity * 3 > 1:
				g = int(min(255 * (intensity * 3 - 1), 255))
			
				if intensity * 3 * 2:
					r = int(min(255 * (intensity * 3 - 2), 255))
		
		colourer.cache[intensity] = (r, g, b)
		
		return(r, g, b)
	
colourer.cache = {}


# These values were sort of determined by experimentation,
# and they may only be suitable for the test region.

# How far it has to stray to be considered escaped
ESCAPE_THRESHOLD = 75

# Maximum number of generations to check for escape.
# Higher number = higher resolution
ESCAPE_TIME_LIMIT = 75

def escape_time(Z):
	"""Returns the escape time (in generations) for a complex point,
	returning None if the point doesn't escape within ESCAPE_THRESHOLD
	generations."""
	
	Z_n = Z
	
	for n in range(ESCAPE_TIME_LIMIT):
		Z_n = (Z_n * Z_n) + Z
		
		if abs(Z - Z_n) > ESCAPE_THRESHOLD:
			break
	
	return(n + 1)

def main(*args):
	dimensions = width, height = 1280 * 2, 800 * 2
	
	x_range, y_range = 4, None
	
	if y_range is None:
		y_range = x_range * height / width
	
	center = -.5, 0
	
	image = Image.new("RGB", dimensions)
	
	pixels = image.load()
	
	data = {}
	
	max_escape_time = 0
	
	start = time.time()
	
	for x in range(width):
		for y in range(height):
			x_adapted = (x / (width - 1) - .5) * x_range + center[0]
			y_adapted = (y / (height - 1) - .5) * y_range + center[1]
			
			this_escape_time = data[x, y] = escape_time(complex(x_adapted, y_adapted))
			
			if this_escape_time > max_escape_time:
				max_escape_time = this_escape_time
			
		completed = (x + 1) / width
		speed = (time.time() - start) / completed
		
		sys.stderr.write("Calculations: {0: >7.2%}, ETA {1:.1f} seconds.\n".format(completed, speed * (1 - completed)))
	sys.stderr.write("Calculations: Completed in {0:.1f} seconds.\n".format(time.time() - start))
	sys.stderr.write("Maximum escape time was {0} generations.\n".format(max_escape_time))
	
	start = time.time()
	
	for x in range(width):
		for y in range(height):
			pixels[x, y] = colourer((data[x, y] - 1) / max_escape_time)
		completed = (x + 1) / width
		
		speed = (time.time() - start) / completed
		sys.stderr.write("Imaging: {0: >7.2%}, ETA {1:.1f} seconds.\n".format(completed, speed * (1 - completed)))
	sys.stderr.write("Imaging: Completed in {0:.1f} seconds.\n".format(time.time() - start))
	
	image.save(os.path.expanduser("~/Desktop/{0}-{1}-{2}x{3}.png".format(ESCAPE_THRESHOLD, ESCAPE_TIME_LIMIT, width, height)))

if __name__ == "__main__": sys.exit(main(*sys.argv[1:]))