|
| 1 | +""" |
| 2 | +=================================== |
| 3 | +Shaded & power normalized rendering |
| 4 | +=================================== |
| 5 | +
|
| 6 | +The Mandelbrot set rendering can be improved by using a normalized recount |
| 7 | +associated with a power normalized colormap (gamma=0.3). Rendering can be |
| 8 | +further enhanced thanks to shading. |
| 9 | +
|
| 10 | +The `maxiter` gives the precision of the computation. `maxiter=200` should |
| 11 | +take a few seconds on most modern laptops. |
| 12 | +""" |
| 13 | +import numpy as np |
| 14 | + |
| 15 | + |
| 16 | +def mandelbrot_set(xmin, xmax, ymin, ymax, xn, yn, maxiter, horizon=2.0): |
| 17 | + X = np.linspace(xmin, xmax, xn, dtype=np.float32) |
| 18 | + Y = np.linspace(ymin, ymax, yn, dtype=np.float32) |
| 19 | + C = X + Y[:, None]*1j |
| 20 | + N = np.zeros(C.shape, dtype=int) |
| 21 | + Z = np.zeros(C.shape, np.complex64) |
| 22 | + for n in range(maxiter): |
| 23 | + I = np.less(abs(Z), horizon) |
| 24 | + N[I] = n |
| 25 | + Z[I] = Z[I]**2 + C[I] |
| 26 | + N[N == maxiter-1] = 0 |
| 27 | + return Z, N |
| 28 | + |
| 29 | + |
| 30 | +if __name__ == '__main__': |
| 31 | + import time |
| 32 | + import matplotlib |
| 33 | + from matplotlib import colors |
| 34 | + import matplotlib.pyplot as plt |
| 35 | + |
| 36 | + xmin, xmax, xn = -2.25, +0.75, 3000/2 |
| 37 | + ymin, ymax, yn = -1.25, +1.25, 2500/2 |
| 38 | + maxiter = 200 |
| 39 | + horizon = 2.0 ** 40 |
| 40 | + log_horizon = np.log(np.log(horizon))/np.log(2) |
| 41 | + Z, N = mandelbrot_set(xmin, xmax, ymin, ymax, xn, yn, maxiter, horizon) |
| 42 | + |
| 43 | + # Normalized recount as explained in: |
| 44 | + # https://linas.org/art-gallery/escape/smooth.html |
| 45 | + # https://www.ibm.com/developerworks/community/blogs/jfp/entry/My_Christmas_Gift |
| 46 | + |
| 47 | + # This line will generate warnings for null values but it is faster to |
| 48 | + # process them afterwards using the nan_to_num |
| 49 | + with np.errstate(invalid='ignore'): |
| 50 | + M = np.nan_to_num(N + 1 - np.log(np.log(abs(Z)))/np.log(2) + log_horizon) |
| 51 | + |
| 52 | + dpi = 72 |
| 53 | + width = 10 |
| 54 | + height = 10*yn/xn |
| 55 | + fig = plt.figure(figsize=(width, height), dpi=dpi) |
| 56 | + ax = fig.add_axes([0.0, 0.0, 1.0, 1.0], frameon=False, aspect=1) |
| 57 | + |
| 58 | + # Shaded rendering |
| 59 | + light = colors.LightSource(azdeg=315, altdeg=10) |
| 60 | + M = light.shade(M, cmap=plt.cm.hot, vert_exag=1.5, |
| 61 | + norm=colors.PowerNorm(0.3), blend_mode='hsv') |
| 62 | + plt.imshow(M, extent=[xmin, xmax, ymin, ymax], interpolation="bicubic") |
| 63 | + ax.set_xticks([]) |
| 64 | + ax.set_yticks([]) |
| 65 | + |
| 66 | + # Some advertisement for matplotlib |
| 67 | + year = time.strftime("%Y") |
| 68 | + major, minor, micro = matplotlib.__version__.split('.', maxsplit=2) |
| 69 | + text = ("The Mandelbrot fractal set\n" |
| 70 | + "Rendered with matplotlib %s.%s, %s — http://matplotlib.org" |
| 71 | + % (major, minor, year)) |
| 72 | + ax.text(xmin+.025, ymin+.025, text, color="white", fontsize=12, alpha=0.5) |
| 73 | + |
| 74 | + plt.show() |
0 commit comments