{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
" \n",
" \n",
" "
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"from lets_plot import *\n",
"\n",
"LetsPlot.setup_html()"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"from numba import jit\n",
"from copy import copy\n",
"from scipy import ndimage"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"@jit\n",
"def log_threshold(z):\n",
" return np.log(np.log(z))/np.log(2)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"@jit\n",
"def mandelbrot(c, threshold=2, num_iter=30, anti_aliasing=True):\n",
" z = c\n",
" for i in range(num_iter):\n",
" if np.abs(z) > threshold:\n",
" return i - (log_threshold(np.abs(z)) - log_threshold(threshold))*anti_aliasing\n",
" z = z*z + c\n",
" return 0"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"@jit\n",
"def mandelbrot_set(x_min, x_max, y_min, y_max, width=270, height=270, threshold=2**40, num_iter=30, power=0.2, anti_aliasing=True):\n",
" xs = np.linspace(x_min, x_max, width)\n",
" ys = np.linspace(y_min, y_max, height)\n",
" img = np.array([mandelbrot(complex(x, y), threshold, num_iter, anti_aliasing)**power for y in ys for x in xs])\n",
" return img"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"@jit\n",
"def partial(Z, axis=0, dx=1.):\n",
" if axis==0:\n",
" return (Z[:-2,1:-1] - Z[2:,1:-1])/dx\n",
" elif axis==1:\n",
" return (Z[1:-1,:-2] - Z[1:-1,2:])/dx"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"@jit\n",
"def laplacian(Z, dx=1.):\n",
" return ( Z[:-2,1:-1] +\n",
" Z[1:-1,:-2] - 4*Z[1:-1,1:-1] + Z[1:-1,2:] +\n",
" Z[2:,1:-1] )/dx**2"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"def get_boundaries(center, span, zoom):\n",
" return center - span/2.**zoom, center + span/2.**zoom "
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"def plot_colored_image(image, rgb, shape, x_size, y_size, title):\n",
" width, height = shape\n",
" colored_image = np.zeros((width, height, 3))\n",
" red, green, blue = rgb\n",
" colored_image[:,:, 0] = image * red\n",
" colored_image[:,:, 1] = image * green\n",
" colored_image[:,:, 2] = image * blue\n",
" plot = ggplot() + geom_image(colored_image) + ggsize(x_size, y_size) + theme(legend_position='none', axis_ticks='blank', axis_line='blank', axis_title='blank', axis_text='blank') + ggtitle(title)\n",
" return plot"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"def plot_mandelbrot_set(zoom=1, center=(-0.5, 0.0), rgb=(0.5, 0.2, 0.9), delta_x=1, delta_y=1, width=500, height=500, threshold=2**40, num_iter_exp=5, power=0.19, x_size=500, y_size=500, dx=1., anti_aliasing=True):\n",
" center_x, center_y = center\n",
" x_min, x_max = get_boundaries(center_x, delta_x, zoom) \n",
" y_min, y_max = get_boundaries(center_y, delta_y, zoom)\n",
" num_iter = 2**num_iter_exp\n",
" image = mandelbrot_set(x_min, x_max, y_min, y_max, width, height, threshold, num_iter, power, anti_aliasing) \n",
" image = image.reshape(width, height)\n",
"\n",
" bunch = GGBunch()\n",
" plot = plot_colored_image(image, rgb, image.shape, x_size, y_size, f'Zoom: {zoom}')\n",
" bunch.add_plot(plot, 0, 0)\n",
"\n",
" diff_image = laplacian(image, dx=dx)\n",
" diff_image = .2*diff_image\n",
" diff_image_shape = (width - 2*int(dx), height - 2*int(dx))\n",
"\n",
" plot = plot_colored_image(diff_image, rgb, diff_image_shape, x_size, y_size, f'Center: {center}')\n",
" bunch.add_plot(plot, x_size, 0)\n",
" bunch.show()"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plot_mandelbrot_set(zoom=16, center=(-0.7487475, 0.0650725), num_iter_exp=12)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Mandelbulb\n",
"There is no strict analogy for complex numbers in 3D space, so there are a number of ways to determine 3D fractals. Some of them use quaternions.\n",
"\n",
"Consider 3D polar coordinates:\n",
"$$\n",
"\\left\\{\n",
"\\begin{array}{l}\n",
"r = \\sqrt{x^2+y^2+z^2},\\\\\n",
"\\theta = \\arctan\\left(\\frac{\\sqrt{x^2+y^2}}{z}\\right),\\\\\n",
"\\phi=\\arctan\\left(\\frac{y}{x}\\right),\n",
"\\end{array}\n",
"\\right.\n",
"$$\n",
"Going back to Descartes coordinates we get:\n",
"$$\\left\\{\n",
"\\begin{array}{l}\n",
"x = r\\sin\\theta\\cos\\phi\\\\\n",
"y = r\\sin\\theta\\sin\\phi\\\\\n",
"z = r\\cos\\theta\n",
"\\end{array}\n",
"\\right.$$\n",
"Given $v = (x,y,z)$, we can define the $n$-th power of this triple:\n",
"$$(x,y,z)^n=r^n\\left(\\sin n\\theta\\cos n\\phi, \\sin n\\theta \\sin n \\phi, \\cos n\\theta\\right)$$\n",
"Now we can define the recurrent formula for the 3D Mandelbrot set (called a Mandelbulb) as $\\displaystyle v_{n+1}=v_n^p + c$, where $c=(x,y,z)$. The typical choice for the power $p$ is 8. \n"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"def get_plane_points(Q, center, span, zoom, width, height):\n",
" x_min, x_max = get_boundaries(center[0], span[0], zoom) \n",
" y_min, y_max = get_boundaries(center[1], span[1], zoom)\n",
" a, b , c = Q\n",
" x = np.linspace(x_min, x_max, width)\n",
" y = np.linspace(y_min, y_max, height)\n",
" x, y = np.meshgrid(x, y)\n",
" x, y = x.reshape(-1), y.reshape(-1)\n",
" if np.abs(c) > 1e-4:\n",
" z = -(a*x + b*y)/c\n",
" P = np.vstack((x, y, z)).T\n",
" elif np.abs(a) > 1e-4:\n",
" z = -(c*x + b*y)/a\n",
" P = np.vstack((z, y, x)).T\n",
" elif np.abs(b) > 1e-4:\n",
" z = -(a*x + c*y)/b\n",
" P = np.vstack((x, z, y)).T\n",
" return P"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"def get_directions(P, Q):\n",
" v = np.array(P - Q)\n",
" v = v/np.linalg.norm(v, axis=1)[:, np.newaxis]\n",
" return v"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"@jit\n",
"def DistanceEstimator(positions, iterations, degree=8, bailout=1000):\n",
" m = positions.shape[0]\n",
" x, y, z = np.zeros(m), np.zeros(m), np.zeros(m)\n",
" x0, y0, z0 = positions[:, 0], positions[:, 1], positions[:, 2]\n",
" dr = np.zeros(m) + 1\n",
" r = np.zeros(m)\n",
" theta = np.zeros(m)\n",
" phi = np.zeros(m)\n",
" zr = np.zeros(m)\n",
" for _ in range(iterations):\n",
" r = np.sqrt(x*x + y*y + z*z)\n",
" idx1 = r < bailout\n",
" dr[idx1] = np.power(r[idx1], degree - 1) * degree * dr[idx1] + 1.0\n",
"\n",
" theta[idx1] = np.arctan2(np.sqrt(x[idx1]*x[idx1] + y[idx1]*y[idx1]), z[idx1])\n",
" phi[idx1] = np.arctan2(y[idx1], x[idx1])\n",
"\n",
" zr[idx1] = r[idx1] ** degree\n",
" theta[idx1] = theta[idx1] * degree\n",
" phi[idx1] = phi[idx1] * degree\n",
"\n",
" x[idx1] = zr[idx1] * np.sin(theta[idx1]) * np.cos(phi[idx1]) + x0[idx1]\n",
" y[idx1] = zr[idx1] * np.sin(theta[idx1]) * np.sin(phi[idx1]) + y0[idx1]\n",
" z[idx1] = zr[idx1] * np.cos(theta[idx1]) + z0[idx1]\n",
"\n",
" return 0.5 * np.log(r) * r / dr"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"def trace(start, directions, max_steps, min_distance, iterations, degree, bailout, power):\n",
" total_distance = np.zeros(directions.shape[0])\n",
" keep_iterations = np.ones_like(total_distance)\n",
" steps = np.zeros_like(total_distance)\n",
" for _ in range(max_steps):\n",
" positions = start[np.newaxis, :] + total_distance[:, np.newaxis] * directions\n",
" distance = DistanceEstimator(positions, iterations, degree, bailout)\n",
" keep_iterations[distance < min_distance] = 0\n",
" total_distance += distance * keep_iterations\n",
" steps += keep_iterations\n",
" return 1 - (steps/max_steps)**power"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"def plot_mandelbulb(degree=8, observer_position=np.array([1., 1., 3.]), max_steps=32, iterations=32, bailout=32000, min_distance=5e-3, zoom=0, power=0.2, width=500, height=500, x_size=500, y_size=500, span=[1.2, 1.2], center=[0, 0]):\n",
" plane_points = get_plane_points(observer_position, center=center, span=span, zoom=zoom, width=width, height=height)\n",
" directions = get_directions(plane_points, observer_position)\n",
" image = trace(observer_position, directions, max_steps, min_distance, iterations, degree, bailout, power)\n",
" image = image.reshape(width, height)\n",
" p = ggplot() \\\n",
" + geom_image(image) + ggsize(x_size, y_size) \\\n",
" + theme(legend_position='none', axis_ticks='blank', axis_line='blank', axis_title='blank', axis_text='blank')\n",
" return p"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"def plot_mandelbulb(degree=8, observer_position=np.array([1., 1., 3.]), max_steps=32, iterations=32, bailout=32000, min_distance=5e-3, zoom=0, power=0.2, width=500, height=500, x_size=500, y_size=500, span=[1.2, 1.2], center=[0, 0]):\n",
" plane_points = get_plane_points(observer_position, center=center, span=span, zoom=zoom, width=width, height=height)\n",
" directions = get_directions(plane_points, observer_position)\n",
" image = trace(observer_position, directions, max_steps, min_distance, iterations, degree, bailout, power)\n",
" image = image.reshape(width, height)\n",
"\n",
" p = ggplot() \\\n",
" + geom_image(image) + ggsize(x_size, y_size) \\\n",
" + theme(legend_position='none', axis_ticks='blank', axis_line='blank', axis_title='blank', axis_text='blank')\n",
"\n",
" return p"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"observer_position = np.array([1., 2., 3.])\n",
"degree = 8 # defines the form of the bulb"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot_mandelbulb(zoom=0, degree=degree, observer_position=observer_position, bailout=2**20, min_distance=5e-4)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot_mandelbulb(zoom=2, degree=degree, observer_position=observer_position, bailout=2**40, min_distance=5e-5)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot_mandelbulb(zoom=3, degree=degree, observer_position=observer_position, bailout=2**40, min_distance=5e-5)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot_mandelbulb(zoom=0, degree=2, observer_position=observer_position, bailout=2**40, min_distance=5e-3)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot_mandelbulb(zoom=0, degree=3, observer_position=observer_position, bailout=2**40, min_distance=5e-3)"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot_mandelbulb(zoom=0, degree=4, observer_position=observer_position, bailout=2**40, min_distance=5e-3)"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot_mandelbulb(zoom=0, degree=5, observer_position=observer_position, bailout=2**40, min_distance=5e-3)"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot_mandelbulb(zoom=0, degree=6, observer_position=observer_position, bailout=2**40, min_distance=5e-3)"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot_mandelbulb(zoom=0, degree=7, observer_position=observer_position, bailout=2**40, min_distance=5e-3)"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot_mandelbulb(zoom=0, degree=8, observer_position=observer_position, bailout=2**40, min_distance=5e-3)"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot_mandelbulb(zoom=0, degree=9, observer_position=observer_position, bailout=2**40, min_distance=5e-3)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}