Mandelbrot set: Difference between revisions

Content added Content deleted
m (→‎Normal Map Effect, Mercator Projection and Deep Zoom Images: Added Jussi Härkönen's paper and made changes to match Arnaud Chéritat's image.)
Line 10,902: Line 10,902:


===Normal Map Effect, Mercator Projection and Deep Zoom Images===
===Normal Map Effect, Mercator Projection and Deep Zoom Images===
The Mandelbrot set is represented by distance estimation and normal maps using NumPy and complex matrices (cf. Arnaud Chéritat: [https://www.math.univ-toulouse.fr/~cheritat/wiki-draw/index.php/Mandelbrot_set#Normal_map_effect ''Normal map effect'']). Note that the second derivative (ddZ) grows very fast, so the second method can only be used for small iteration numbers (n <= 400). See also [https://www.shadertoy.com/view/wtscDX Julia Stripes] on Shadertoy.
The Mandelbrot set is represented by distance estimation and normal maps using NumPy and complex matrices (cf. Arnaud Chéritat: [https://www.math.univ-toulouse.fr/~cheritat/wiki-draw/index.php/Mandelbrot_set#Normal_map_effect ''Normal map effect'']). Note that the second derivative (ddZ) grows very fast, so the second method can only be used for small iteration numbers (n <= 400). See also [https://web.archive.org/web/20140618000747/http://jussiharkonen.com/files/on_fractal_coloring_techniques(lo-res).pdf ''On Smooth Fractal Coloring Techniques''] by Jussi Härkönen and [https://www.shadertoy.com/view/wtscDX Julia Stripes] on Shadertoy.
<syntaxhighlight lang="python">import numpy as np
<syntaxhighlight lang="python">import numpy as np
import matplotlib.pyplot as plt
import matplotlib.pyplot as plt
Line 10,910: Line 10,910:


direction, height = 45, 1.5 # direction and height of the incoming light
direction, height = 45, 1.5 # direction and height of the incoming light
stripes, damping = 10, 2.0 # stripe density and damping parameter
stripes, damping = 4, 2.0 # stripe density and damping parameter


x = np.linspace(0, 2, num=d+1)
x = np.linspace(0, 2, num=d+1)
Line 10,922: Line 10,922:


for k in range(n):
for k in range(n):
M = Z.real ** 2 + Z.imag ** 2 < r ** 2
M = abs(Z) < r
S[M], T[M] = S[M] + np.sin(stripes * np.angle(Z[M])), T[M] + 1
S[M], T[M] = S[M] + np.cos(stripes * np.angle(Z[M])), T[M] + 1
Z[M], dZ[M], ddZ[M] = Z[M] ** 2 + C[M], 2 * Z[M] * dZ[M] + 1, 2 * (dZ[M] ** 2 + Z[M] * ddZ[M])
Z[M], dZ[M], ddZ[M] = Z[M] ** 2 + C[M], 2 * Z[M] * dZ[M] + 1, 2 * (dZ[M] ** 2 + Z[M] * ddZ[M])


N = abs(Z) >= r # normal map effect 1 (potential function)
N = abs(Z) >= r # normal map effect 1 (potential function)
P, Q = S[N] / T[N], (S[N] + np.sin(stripes * np.angle(Z[N]))) / (T[N] + 1)
P, Q = S[N] / T[N], (S[N] + np.cos(stripes * np.angle(Z[N]))) / (T[N] + 1)
F = np.log2(np.log(np.abs(Z[N])) / np.log(r)) # fraction between 0 and 1 (for interpolation)
F = np.log2(np.log(np.abs(Z[N])) / np.log(r)) # fraction between 0 and 1 (for interpolation)
H = F * P + (1 - F) * Q # height perturbation (by linear interpolation)
R = Q + (P - Q) * F * F * (3 - 2 * F) # hermite interpolation (r is between q and p)
U = Z[N] / dZ[N] # normal vectors to the equipotential lines
U, H = Z[N] / dZ[N], 1 + R / damping # normal vectors to the equipotential lines and height perturbation
U, v = U / abs(U), np.exp(direction / 180 * np.pi * 1j) # unit normal vectors and unit 2D vector
U, v = U / abs(U), np.exp(direction / 180 * np.pi * 1j) # unit normal vectors and vector in light direction
D[N] = np.maximum((U.real * v.real + U.imag * v.imag + (1 + H / damping) * height) / (1 + height), 0)
D[N] = np.maximum((U.real * v.real + U.imag * v.imag + H * height) / (1 + height), 0)


plt.imshow(D ** 1.0, cmap=plt.cm.bone, origin="lower")
plt.imshow(D ** 1.0, cmap=plt.cm.bone, origin="lower")
Line 10,939: Line 10,939:
N = abs(Z) >= r # normal map effect 2 (distance estimation)
N = abs(Z) >= r # normal map effect 2 (distance estimation)
U = Z[N] * dZ[N] * ((1 + np.log(abs(Z[N]))) * np.conj(dZ[N] ** 2) - np.log(abs(Z[N])) * np.conj(Z[N] * ddZ[N]))
U = Z[N] * dZ[N] * ((1 + np.log(abs(Z[N]))) * np.conj(dZ[N] ** 2) - np.log(abs(Z[N])) * np.conj(Z[N] * ddZ[N]))
U, v = U / abs(U), np.exp(direction / 180 * np.pi * 1j) # unit normal vectors and unit 2D vector
U, v = U / abs(U), np.exp(direction / 180 * np.pi * 1j) # unit normal vectors and vector in light direction
D[N] = np.maximum((U.real * v.real + U.imag * v.imag + height) / (1 + height), 0)
D[N] = np.maximum((U.real * v.real + U.imag * v.imag + height) / (1 + height), 0)