{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "
\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": [ "