代碼：46 lines (34 sloc) 1.01 kb
'''
a fast mandelbrot set wallpaper renderer
reddit discussion: https://www.reddit.com/r/math/comments/2abwyt/smooth_colour_mandelbrot/
import numpy as np
from pil import image
from numba import jit
maxiters = 200
radius = 100
@jit
def color(z, i):
v = np.log2(i + 1 - np.log2(np.log2(abs(z)))) / 5
if v < 1.0:
return v**4, v**2.5, v
else:
v = max(0, 2-v)
return v, v**1.5, v**3
def iterate(c):
z = 0j
for i in range(maxiters):
if z.real*z.real + z.imag*z.imag > radius:
return color(z, i)
z = z*z + c
return 0, 0 ,0
def main(xmin, xmax, ymin, ymax, width, height):
x = np.linspace(xmin, xmax, width)
y = np.linspace(ymax, ymin, height)
z = x[none, :] + y[:, none]*1j
red, green, blue = np.asarray(np.frompyfunc(iterate, 1, 3)(z)).astype(np.float)
img = np.dstack((red, green, blue))
image.fromarray(np.uint8(img*255)).save('mandelbrot.png')
if __name__ == '__main__':
main(-2.1, 0.8, -1.16, 1.16, 1200, 960)
多米諾洗牌算法
代碼連結：https://github.com/neozhaoliang/pywonderland/tree/master/src/domino
正二十面體萬花筒
代碼：53 lines (40 sloc) 1.24 kb
a kaleidoscope pattern with icosahedral symmetry.
from matplotlib.colors import hsv_to_rgb
def klein(z):
'''klein's j-function'''
return 1728 * (z * (z**10 + 11 * z**5 - 1))**5 / \
(-(z**20 + 1) + 228 * (z**15 - z**5) - 494 * z**10)**3
def riemannsphere(z):
map the complex plane to riemann's sphere via stereographic projection
'''
t = 1 + z.real*z.real + z.imag*z.imag
return 2*z.real/t, 2*z.imag/t, 2/t-1
def mobius(z):
distort the result image by a mobius transformation
return (z - 20)/(3*z + 1j)
def main(imgsize):
x = np.linspace(-6, 6, imgsize)
y = np.linspace(6, -6, imgsize)
z = riemannsphere(klein(mobius(klein(z))))
# define colors in hsv space
h = np.sin(z[0]*np.pi)**2
s = np.cos(z[1]*np.pi)**2
v = abs(np.sin(z[2]*np.pi) * np.cos(z[2]*np.pi))**0.2
hsv = np.dstack((h, s, v))
# transform to rgb space
img = hsv_to_rgb(hsv)
image.fromarray(np.uint8(img*255)).save('kaleidoscope.png')
import time
start = time.time()
main(imgsize=800)
end = time.time()
print('runtime: {:3f} seconds'.format(end - start))
newton 疊代分形
代碼：46 lines (35 sloc) 1.05 kb
import matplotlib.pyplot as plt
# define functions manually, do not use numpy's poly1d funciton!
@jit('complex64(complex64)', nopython=true)
def f(z):
# z*z*z is faster than z**3
return z*z*z - 1
def df(z):
return 3*z*z
@jit('float64(complex64)', nopython=true)
def iterate(z):
num = 0
while abs(f(z)) > 1e-4:
w = z - f(z)/df(z)
num += np.exp(-1/abs(w-z))
z = w
return num
def render(imgsize):
x = np.linspace(-1, 1, imgsize)
y = np.linspace(1, -1, imgsize)
z = x[none, :] + y[:, none] * 1j
img = np.frompyfunc(iterate, 1, 1)(z).astype(np.float)
fig = plt.figure(figsize=(imgsize/100.0, imgsize/100.0), dpi=100)
ax = fig.add_axes([0, 0, 1, 1], aspect=1)
ax.axis('off')
ax.imshow(img, cmap='hot')
fig.savefig('newton.png')
render(imgsize=400)
print('runtime: {:03f} seconds'.format(end - start))
李代數e8 的根系
代碼連結：https://github.com/neozhaoliang/pywonderland/blob/master/src/misc/e8.py
模群的基本域
代碼連結：
https://github.com/neozhaoliang/pywonderland/blob/master/src/misc/modulargroup.py
彭羅斯鋪砌
https://github.com/neozhaoliang/pywonderland/blob/master/src/misc/penrose.py
wilson 算法
代碼連結：https://github.com/neozhaoliang/pywonderland/tree/master/src/wilson
反應擴散方程模拟
代碼連結：https://github.com/neozhaoliang/pywonderland/tree/master/src/grayscott
120 胞腔
代碼：69 lines (48 sloc) 2.18 kb
# pylint: disable=unused-import
# pylint: disable=undefined-variable
from itertools import combinations, product
from vapory import *
class penrose(object):
grids = [np.exp(2j * np.pi * i / 5) for i in range(5)]
def __init__(self, num_lines, shift, thin_color, fat_color, **config):
self.num_lines = num_lines
self.shift = shift
self.thin_color = thin_color
self.fat_color = fat_color
self.objs = self.compute_pov_objs(**config)
def compute_pov_objs(self, **config):
objects_pool = []
for rhombi, color in self.tile():
p1, p2, p3, p4 = rhombi
polygon = polygon(5, p1, p2, p3, p4, p1,
texture(pigment('color', color), config['default']))
objects_pool.append(polygon)
for p, q in zip(rhombi, [p2, p3, p4, p1]):
cylinder = cylinder(p, q, config['edge_thickness'], config['edge_texture'])
objects_pool.append(cylinder)
for point in rhombi:
x, y = point
sphere = sphere((x, y, 0), config['vertex_size'], config['vertex_texture'])
objects_pool.append(sphere)
return object(union(*objects_pool))
def rhombus(self, r, s, kr, ks):
if (s - r)**2 % 5 == 1:
color = self.thin_color
color = self.fat_color
point = (penrose.grids[r] * (ks - self.shift[s])
- penrose.grids[s] * (kr - self.shift[r])) *1j / penrose.grids[s-r].imag
index = [np.ceil((point/grid).real + shift)
for grid, shift in zip(penrose.grids, self.shift)]
vertices = []
for index[r], index[s] in [(kr, ks), (kr+1, ks), (kr+1, ks+1), (kr, ks+1)]:
vertices.append(np.dot(index, penrose.grids))
vertices_real = [(z.real, z.imag) for z in vertices]
return vertices_real, color
def tile(self):
for r, s in combinations(range(5), 2):
for kr, ks in product(range(-self.num_lines, self.num_lines+1), repeat=2):
yield self.rhombus(r, s, kr, ks)
def put_objs(self, *args):
return object(self.objs, *args)
原文釋出時間為：2017-04-15
本文來自雲栖社群合作夥伴“大資料文摘”，了解相關資訊可以關注“bigdatadigest”微信公衆号