-
Notifications
You must be signed in to change notification settings - Fork 0
/
fractal-images.py
199 lines (193 loc) · 7.44 KB
/
fractal-images.py
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
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
from numpy import *
def mandel(n, m, itermax, xmin, xmax, ymin, ymax):
'''
Fast mandelbrot computation using numpy.
(n, m) are the output image dimensions
itermax is the maximum number of iterations to do
xmin, xmax, ymin, ymax specify the region of the
set to compute.
'''
# The point of ix and iy is that they are 2D arrays
# giving the x-coord and y-coord at each point in
# the array. The reason for doing this will become
# clear below...
ix, iy = mgrid[0:n, 0:m]
# Now x and y are the x-values and y-values at each
# point in the array, linspace(start, end, n)
# is an array of n linearly spaced points between
# start and end, and we then index this array using
# numpy fancy indexing. If A is an array and I is
# an array of indices, then A[I] has the same shape
# as I and at each place i in I has the value A[i].
x = linspace(xmin, xmax, n)[ix]
y = linspace(ymin, ymax, m)[iy]
# c is the complex number with the given x, y coords
c = x+complex(0,1)*y
del x, y # save a bit of memory, we only need z
# the output image coloured according to the number
# of iterations it takes to get to the boundary
# abs(z)>2
img = zeros(c.shape, dtype=int)
# Here is where the improvement over the standard
# algorithm for drawing fractals in numpy comes in.
# We flatten all the arrays ix, iy and c. This
# flattening doesn't use any more memory because
# we are just changing the shape of the array, the
# data in memory stays the same. It also affects
# each array in the same way, so that index i in
# array c has x, y coords ix[i], iy[i]. The way the
# algorithm works is that whenever abs(z)>2 we
# remove the corresponding index from each of the
# arrays ix, iy and c. Since we do the same thing
# to each array, the correspondence between c and
# the x, y coords stored in ix and iy is kept.
ix.shape = n*m
iy.shape = n*m
c.shape = n*m
# we iterate z->z^2+c with z starting at 0, but the
# first iteration makes z=c so we just start there.
# We need to copy c because otherwise the operation
# z->z^2 will send c->c^2.
z = copy(c)
for i in xrange(itermax):
if not len(z): break # all points have escaped
# equivalent to z = z*z+c but quicker and uses
# less memory
multiply(z, z, z)
add(z, c, z)
# these are the points that have escaped
rem = abs(z)>2.0
# colour them with the iteration number, we
# add one so that points which haven't
# escaped have 0 as their iteration number,
# this is why we keep the arrays ix and iy
# because we need to know which point in img
# to colour
img[ix[rem], iy[rem]] = i+1
# -rem is the array of points which haven't
# escaped, in numpy -A for a boolean array A
# is the NOT operation.
rem = -rem
# So we select out the points in
# z, ix, iy and c which are still to be
# iterated on in the next step
z = z[rem]
ix, iy = ix[rem], iy[rem]
c = c[rem]
return img
def julia(n, m, itermax, xmin, xmax, ymin, ymax, name):
'''
Fast mandelbrot computation using numpy.
(n, m) are the output image dimensions
itermax is the maximum number of iterations to do
xmin, xmax, ymin, ymax specify the region of the
set to compute.
'''
# The point of ix and iy is that they are 2D arrays
# giving the x-coord and y-coord at each point in
# the array. The reason for doing this will become
# clear below...
ix, iy = mgrid[0:n, 0:m]
# Now x and y are the x-values and y-values at each
# point in the array, linspace(start, end, n)
# is an array of n linearly spaced points between
# start and end, and we then index this array using
# numpy fancy indexing. If A is an array and I is
# an array of indices, then A[I] has the same shape
# as I and at each place i in I has the value A[i].
x = linspace(xmin, xmax, n)[ix]
y = linspace(ymin, ymax, m)[iy]
# z is the complex number with the given x, y coords
z = x + complex(0,1)*y
del x, y # save a bit of memory, we only need z
# c = random(1)*1.2 - .5 + complex(0,1)*(random(1)*1.-.5)
r = 0.8 + random(1) / 5
phi = random(1)*2*pi
c = r*cos(phi) + complex(0,1)*r*sin(phi)
# the output image coloured according to the number
# of iterations it takes to get to the boundary
# abs(z)>2
img = zeros(z.shape, dtype=int)
# Here is where the improvement over the standard
# algorithm for drawing fractals in numpy comes in.
# We flatten all the arrays ix, iy and c. This
# flattening doesn't use any more memory because
# we are just changing the shape of the array, the
# data in memory stays the same. It also affects
# each array in the same way, so that index i in
# array c has x, y coords ix[i], iy[i]. The way the
# algorithm works is that whenever abs(z)>2 we
# remove the corresponding index from each of the
# arrays ix, iy and c. Since we do the same thing
# to each array, the correspondence between c and
# the x, y coords stored in ix and iy is kept.
ix.shape = n*m
iy.shape = n*m
z.shape = n*m
# we iterate z->z^2+c with z starting at 0, but the
# first iteration makes z=c so we just start there.
# We need to copy c because otherwise the operation
# z->z^2 will send c->c^2.
# z = copy(c)
PARS = random(3) / 2
for i in xrange(itermax):
if not len(z): break # all points have escaped
# equivalent to z = z*z+c but quicker and uses
# less memory
if random(1) < PARS[0] :
z = log(z)
if random(1) < PARS[1] :
multiply(z, z, z)
if random(1) < 0.5 :
z = 1 / (complex(1,1)-z)
else :
multiply(z, z, z)
if random(1) < PARS[2] :
add(z, c, z)
else :
add(z, -c, z)
# these are the points that have escaped
rem = abs(z)>2.0
# colour them with the iteration number, we
# add one so that points which haven't
# escaped have 0 as their iteration number,
# this is why we keep the arrays ix and iy
# because we need to know which point in img
# to colour
img[ix[rem], iy[rem]] = i+1
# -rem is the array of points which haven't
# escaped, in numpy -A for a boolean array A
# is the NOT operation.
rem = ~rem
# So we select out the points in
# z, ix, iy and c which are still to be
# iterated on in the next step
z = z[rem]
ix, iy = ix[rem], iy[rem]
# c = c[rem]
timg = log(img+1)
#timg[timg==0] = 101
if( i > 10 ):
continue
ttimg = imshow(timg.T, origin='lower left')
ttimg.write_png(name+'_'+str(i)+'.png')
# ttimg.write_png(name+'_'+str(i)+'.png', noscale=True)
return timg
if __name__=='__main__':
from pylab import *
import time
import sys
if( len(sys.argv) > 1 ):
name = str(sys.argv[1])
else:
name = 'fc'
for ITER in range(0, 1):
# print ITER
# start = time.time()
I = julia(800, 800, 100, -2, 2, -2, 2, name)
# print 'Time taken:', time.time()-start
I[I==0] = 101
img = imshow(I.T, origin='lower left')
img.write_png(name+'.png')
del I, img
#show()