Generalized Apollonian Gasket Fractal (Python recipe) by FB36

ActiveState Code (http://code.activestate.com/recipes/578029/)

# Generalized Apollonian Gasket Fractal
# http://en.wikipedia.org/wiki/Apollonian_Gasket
# http://en.wikipedia.org/wiki/Descartes'_theorem
# http://en.wikibooks.org/wiki/Fractals/Apollonian_fractals
# FB - 20120131
import math
import cmath
import random
from PIL import Image, ImageDraw
imgx = 512
imgy = 512
image = Image.new("RGB", (imgx, imgy))
draw = ImageDraw.Draw(image)
maxIt = 2000 # of iterations
n = random.randint(2, 7) # of main circles > 1
print 'Number of main circles: ' + str(n)

eps = 0.0001
def appEq(u, v): # approximately equal excluding signs
    return abs(abs(u) - abs(v)) <= eps

# test if 2 circles are approximately equal
def appEqCir(a, b):
    return (abs(a[0] - b[0]) <= eps and appEq(a[1], b[1]))

# test if 2 circles are tangent
def isTanCir(a, b):
    return (appEq(a[0] - b[0], abs(a[1]) + abs(b[1])) or appEq(a[0] - b[0], abs(a[1]) - abs(b[1])))

# test if 2 circles are intersecting
def isIntCir(a, b):
    dist = abs(a[0] - b[0]) # distance between 2 centers
    rmin = min(abs(a[1]), abs(b[1])) # min radius
    rmax = max(abs(a[1]), abs(b[1])) # max radius
    return (dist + eps < rmin + rmax and dist + rmin > rmax + eps)

# test if 3 circles are tangent to each other
def isTangent(a, b, c):
    flag = True
    if appEqCir(a, b) or appEqCir(a, c) or appEqCir(b, c): # all circles must be different
        flag = False
    if not isTanCir(a, b):
        flag = False
    if not isTanCir(a, c):
        flag = False
    if not isTanCir(b, c):
        flag = False
    return flag

# test if the 4th circle is tangent to previous 3
def isTan4(a, b, c, d):
    return isTangent(a, b, d) and isTangent(a, c, d) and isTangent(b, c, d)

# input: 3 circles w/ each as a tuple: (complex(x,y),r)
# output: 4 tangent circles as a tuple
def cif(a, b, c):
    # curvatures of the circles
    k1 = 1.0 / a[1]
    k2 = 1.0 / b[1]
    k3 = 1.0 / c[1]
    # curvatures of tangent circles
    temp0 = 2.0 * math.sqrt(k1 * k2 + k1 * k3 + k2 * k3)
    temp1 = k1 + k2 + k3
    k40 = temp1 + temp0
    k41 = temp1 - temp0
    # centers of tangent circles
    temp0 = 2.0 * cmath.sqrt(k1*k2*a[0]*b[0]+k1*k3*a[0]*c[0]+k2*k3*b[0]*c[0])
    temp1 = k1 * a[0] + k2 * b[0] + k3 * c[0]
    c40 = (temp1 + temp0) / k40
    c41 = (temp1 - temp0) / k41
    c42 = (temp1 - temp0) / k40
    c43 = (temp1 + temp0) / k41
    # radius of tangent circles
    r0 = 1.0 / k40
    r1 = 1.0 / k41
    # there are 4 solutions and only 2 are correct (tangent) circles
    # those are c0 and c1 or c2 and c3
    c0 = (c40, r0)
    c1 = (c41, r1)
    c2 = (c42, r0)
    c3 = (c43, r1)
    if isTan4(a, b, c, c0) and isTan4(a, b, c, c1):
        return (c0, c1)
    else:
        return (c2, c3)

# pixel coordinate calculation
def pxy(x, y):
    kx = int((x - xa) / (xb - xa) * (imgx - 1))
    ky = int((y - ya) / (yb - ya) * (imgy - 1))
    return (kx, ky)   

# generate the main circles: n circles in a circle
circles = []
a = 2.0 * math.pi / n
# radius of the inner circles
r = math.hypot(math.cos(a) - math.cos(0), math.sin(a) - math.sin(0)) / 2.0
r0 = 1.0 + r # radius of the outer circle
r1 = r0 - 2.0 * r # radius of the center circle
# origin
ox = r0
oy = r0
rmin = 2.0 * r0 / math.hypot(imgx, imgy) # min radius that can be drawn
# drawing area
xa = ox - r0
xb = ox + r0
ya = oy - r0
yb = oy + r0

circles.append((complex(ox, oy), -r0)) # add the outer circle w/ - curvature
for j in range(n): # add the main circles
    cx = math.cos(a * j) + ox
    cy = math.sin(a * j) + oy
    circles.append((complex(cx, cy), r))
if r1 > rmin: # if possible add a center circle
    circles.append((complex(ox, oy), r1))    

p0 = 0
for i in range(maxIt):
    p = int(100 * i / (maxIt - 1))
    if p > p0:
        p0 = p
        print "%" + str(p)
    # randomly select 3 circles
    while True:
        a = random.choice(circles)
        b = random.choice(circles)
        c = random.choice(circles)
        # test if all 3 circles are tangent to each other
        if isTangent(a, b, c): 
            break
    # find inversion circles
    d = cif(a, b, c)
    for e in d:
        dx = e[0].real
        dy = e[0].imag
        dr = e[1]
        # ignore the outer circle (- curvature) solutions and too small circles
        if dx > 0.0 and dy > 0.0 and dr > rmin:
            # add the new circle if it was not added before and not intersecting
            flag = True
            for cir in circles:
                # new circle must not be equal/intersecting to any existing circle
                if appEqCir(cir, e) or isIntCir(cir, e):
                    flag = False
                    break
            if flag:
                circles.append(e)

# draw all circles
for cir in circles:
    dx = cir[0].real
    dy = cir[0].imag
    dr = abs(cir[1])
    try:
        draw.ellipse((pxy(dx - dr, dy - dr), pxy(dx + dr, dy + dr)))
    except:
        pass
image.save("ApollonianGasket_" + str(n) + ".png", "PNG")