[toc]
初衷
曼德布洛特集合是在李善友老師的混沌大學(xué)公開(kāi)課上聽(tīng)到的，在極客學(xué)院使用TensorFlow也可以實(shí)現(xiàn)該集合炎咖，將其可視化時(shí)赃泡，于是嘗試模擬該模型。
參考文檔
Numpy乘盼、matplotlib實(shí)現(xiàn)二維數(shù)據(jù)到圖像的轉(zhuǎn)換升熊，添加colormap，無(wú)邊距顯示：https://blog.csdn.net/u010105243/article/details/76695400
極客學(xué)院源碼曼德布洛特集合源碼（較老）：http://wiki.jikexueyuan.com/project/tensorflow-zh/tutorials/mandelbrot.html
源碼實(shí)現(xiàn)曼德勃羅集合
"""
本案例來(lái)展示偏微分方程
使用tensorflow來(lái)實(shí)現(xiàn)該模型
"""
import tensorflow as tf
import numpy as np
import PIL.Image
import PIL.ImageDraw
import matplotlib.pyplot as plt
def DisplayFractal(a, fmt='jpeg'):
"""顯示迭代計(jì)算出的彩色分形圖像绸栅。"""
a_cyclic = (6.28*a/20.0).reshape(list(a.shape)+[1])
img = np.concatenate([10+20*np.cos(a_cyclic),
30+50*np.sin(a_cyclic),
155-80*np.cos(a_cyclic)], 2)
img[a==a.max()] = 0
a = img
a = np.uint8(np.clip(a, 0, 255))
img1 = PIL.Image.fromarray(a)
plt.imsave("image_tf.png", img1)
plt.show()
sess = tf.InteractiveSession()
# 使用NumPy創(chuàng)建一個(gè)在[-2,2]x[-2,2]范圍內(nèi)的2維復(fù)數(shù)數(shù)組
Y, X = np.mgrid[-1.3:1.3:0.005, -2:1:0.005]
Z = X+1j*Y
xs = tf.constant(Z.astype("complex64"))
zs = tf.Variable(xs)
ns = tf.Variable(tf.zeros_like(xs, "float32"))
tf.global_variables_initializer().run()
# 計(jì)算一個(gè)新值z(mì): z^2 + x
zs_ = zs*zs + xs
# 這個(gè)新值會(huì)發(fā)散嗎僚碎？
not_diverged = tf.abs(zs_) < 4
# 更新zs并且迭代計(jì)算。
#
# 說(shuō)明：在這些值發(fā)散之后阴幌，我們?nèi)匀辉谟?jì)算zs，這個(gè)計(jì)算消耗特別大卷中！
# 如果稍微簡(jiǎn)單點(diǎn)矛双，這里有更好的方法來(lái)處理。
#
step = tf.group(
zs.assign(zs_),
ns.assign_add(tf.cast(not_diverged, "float32"))
)
for i in range(200):
step.run()
DisplayFractal(ns.eval())
這里對(duì)wiki上的源碼做了三處改動(dòng)
1. import matplotlib.pyplot as plt
使用plt繪制二維數(shù)組圖
2. img1 = PIL.Image.fromarray(a)
plt.imsave("image_tf.png", img1)
用plt繪圖蟆豫，而不是IPython.display
3. # 這個(gè)新值會(huì)發(fā)散嗎议忽？
not_diverged = tf.abs(zs_) < 4
原函數(shù)tf.complex_abs已經(jīng)在新版tf中棄用。
結(jié)果如下十减，在當(dāng)前目錄下生成一張圖片image_tf.png：
image_tf.png
后記
還有其他繪制曼德勃羅集合的方式栈幸，這里就不展開(kāi)了愤估。
"""
本案例用來(lái)展示曼德布洛特迭代模型.
使用純粹的python 算法實(shí)現(xiàn)該模型, machine_10 使用tensorflow來(lái)實(shí)現(xiàn)該模型
繪圖使用PIL
"""
import time
from PIL import Image, ImageDraw
g_size = (400, 300) # 圖形最終尺寸
g_max_iteration = 256 # 最大迭代次數(shù)
g_bailout = 4 # 最大域
g_zoom = 2.5 / g_size[0] # 縮放參數(shù)
g_offset = (-g_size[0] * 0.25, 0) # 偏移量
g_HSL = (210, 80, 50) # HSL色彩基調(diào)
def draw(antialias = True):
zi = 2 if antialias else 1 # antialias: 抗鋸齒 size = [i * zi
size = [i * zi for i in g_size]
zoom = g_zoom / zi
offset = [i * zi for i in g_offset]
bailout = g_bailout * zi
img = Image.new("RGB", size, 0xffffff)
dr = ImageDraw.Draw(img)
print("painting Mandelbrot Set..")
for xy, color in getPoints(size, offset, zoom):
dr.point(xy, fill = color)
print("100%n")
del dr
if antialias:
img = img.resize(g_size, Image.ANTIALIAS)
img.show()
img.save("mandelbrot_set_%dx%d.png" % g_size)
def getPoints(size, offset, zoom, ti = 0, tstep = 1):
"生成需要繪制的點(diǎn)的坐標(biāo)及顏色"
def getRepeats(c):
z = c
repeats = 0
while abs(z) < g_bailout and repeats < g_max_iteration:
z = z * z + c
repeats += 1
return repeats
def getColor(r):
color = "hsl(0, 0%, 0%)"
if r < g_max_iteration:
v = 1.0 * r / g_max_iteration
h = ch * (1 - v)
s = cs
l = cl * (1 + v)
color = "hsl(%d, %d%%, %d%%)" % (h, s, l)
return color
xs, ys = size
xw, yh = xs / 2, ys / 2
xo, yo = offset
ch, cs, cl = g_HSL
progress = 0
for iy in range(ys):
p = iy * 100 / ys
if iy % 10 == 0 and p != progress:
print ("%d%%..." % p) # 顯示進(jìn)度
progress = p
for ix in range(ti, xs, tstep):
x = (ix - xw + xo) * zoom
y = (iy - yh + yo) * zoom
c = complex(x, y)
r = getRepeats(c)
yield (ix, iy), getColor(r)
def main():
t0 = time.time()
draw()
t = time.time() - t0
print("%dm%.3fs" % (t / 60, t % 60))
if __name__ == "__main__":
main()
