Welcome, guest | Sign In | My Account | Store | Cart

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.

Python, 159 lines
  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")

9 comments

yishu 5 years, 1 month ago  # | flag

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!!

FB36 (author) 5 years, 1 month ago  # | flag

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.

yishu 5 years ago  # | flag

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!

FB36 (author) 5 years ago  # | flag

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.

yishu 5 years ago  # | flag

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?

FB36 (author) 5 years ago  # | flag

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

yishu 5 years ago  # | flag

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.

FB36 (author) 5 years ago  # | flag

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

John Park 4 years, 11 months ago  # | flag

There are some missing circles after generation..