It draws Apollonian Gasket Fractal for any n using Descartes theorem.

This is not the standard way though. It simply randomly finds 3 tangent circles at each iteration and tries to add new circles. The good thing is it can start w/ any arbitrary configuration of main circles, unlike the standard way. The bad thing is some circles will stay missing because of random selections they never get a chance. You can increase the maxIt to get better result but it would slow it down a lot.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 | ```
# 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")
``` |

thank you so much for your code! However, there is one equation I have not understand. Could you explain the equation? How to understand this equation? r = math.hypot(math.cos(a) - math.cos(0), math.sin(a) - math.sin(0)) / 2.0

Thank you so much!!

If you put N circles (as big as possible) of same size inside a bigger circle, that equation gives you the radius of each smaller circles. It is the distance between 2 centers of circles inside (next to each other) divided by 2.

Thank you for your answer. But I still cannot figure out the geometry realtion here. Do you mind to explain it more clearly?

Besides, in your code, the generation of each circle is random. Do you have any idea on generating circles by generations, which means generating circles from big to small?

Thanks again!

I don't think I can explain better without drawings etc also it was years ago for me. But I think my idea was instead of starting with a unit circle and trying to figure out how to pack N circles inside of it, I realized the problem would be simpler if the N circles all have centers on the circumference of unit circle (not inside of it). That equation comes from that. You need to a drawing to see it. And once you have N circles, it it not hard to calculate the radius of a circle that surrounds those N circles and the radius of the small circle in the center of N circles.

Generating the circles in generations/levels is something I couldn't figure out myself. That is why I made it random. But generating them randomly of course is not a good efficient way. That is why I gave up on this method and switched to a much better way later. I don't know if you tried my other programs in this website titled "Circle Inversion Fractals". They also draw Apollonian Gasket fractals but in a much better way.

Thank you for your answer. I got it now!

I found a paper that introduced the packing method generation by generation. http://www.math.columbia.edu/~staff/RTGPapers/ApollonianSums.pdf

It is not easy for me to write my own packing code because I am a beginner just startin learning the packing method. I noticed you assign colors (or images) to different circles, is there any special tips here?

I think that paper you found has the solutions for both creating generations and coloring. (But it requires spending sometime to understand :-)

I usually avoid using any advanced coloring methods for any fractals to keep the code short and easy to understand. So I use either randomly assigned colors or just create colors using mod operation:

if k is fractal iteration number/level etc then:

red=(k%a) * (256/a)

green=(k%b) * (256/b)

blue=(k%c) * (256/c)

a,b,c values selected from 4,8,16,32,64,128 and a<>b,a<>c,b<>c

Thank you so much for your patience! I will take time to read the papers and try to understrand the coloring methods.

In 2D packing ,the ceter is straightforward and can be thought of as a point in the 2D-plane and be described by a complex number. Do you have any idea on how to represent the center of sphere in 3D packing? I am trying to find a way, but cannot figure it out now.

Since there are no 3d extension of complex numbers, I am guessing 3d solution must use 4d numbers: https://en.wikipedia.org/wiki/Quaternion

There are some missing circles after generation..