23.10.2025 13:58
Andrey-82 7 2100 Источник

Привет, Хабр!

Меня всё также зовут Андрей Гринблат. В прошлых материалах я рассказывал о построении фотореалистичных изображений трёхмерных фракталов (часть 1 и часть 2). Это — завершающая статья цикла, в ней я разберу визуализацию оболочки Мандельброта, четырёхмерных аналогов множеств Мандельброта и Жюлиа, и рассмотрю гибридные фракталы.

Оболочка Мандельброта (Mandelbulb)

В 2009 году Даниел Уайт и Пол Ниландер представили трёхмерный аналог множества Мандельброта — Mandelbulb. Их идея состояла в том, чтобы использовать формулу          z_{n+1} = z_n^k+c не для комплексных чисел, а для «трёхмерных».

«Хороших» триплексных чисел нет. Но можно использовать сферические координаты и по аналогии с комплексными числами определить операцию умножения и возведения в степень точек трёхмерного пространства.

\begin{aligned}&p = (x, y, z) \;\longrightarrow\; (r, \theta, \phi), \\[6pt]&\text{где:} \\[4pt]&r = \sqrt{x^2 + y^2 + z^2},\qquad \theta = \arccos\!\left(\dfrac{z}{r}\right), \qquad \phi = \arctan\!\left(\dfrac{y}{x}\right), \\[8pt]&\text{тогда определим возведение в степень:} \\[4pt]&p^k = (x, y, z)^k \;\longrightarrow\; (r, \theta, \phi)^k = (r^k,\; k\theta,\; k\phi), \\[8pt]&\text{переходя обратно в декартовы координаты, получаем:} \\[4pt]&p^k = r^k \!\cdot\!\bigl(\sin(k\theta)\cos(k\phi),\;\sin(k\theta)\sin(k\phi),\;\cos(k\theta)\bigr).\end{aligned}

И теперь скажем (по определению), что точка p принадлежит фракталу, если её орбита, определяемая уравнением p_{n+1} = p_n^k + c не уходит на бесконечность.

Давайте теперь найдём расстояние до фрактала.

Для фракталов, определяемых через «время убегания», хорошо подходит оценка расстояния с помощью якобиана, используемого преобразования, или оценка через потенциал Грина. В нашем случае, аналитически вычислять якобиан не просто. Более удобной и быстро вычисляемой будет оценка расстояния, полученная с помощью логарифмической потенциал-функции (например, аналогичной потенциалу Хаббарда-Дуади). Идея потенциала состоит в том, что это функция, ведущая себя подобно функции расстояния, особенно вблизи фрактала, но в отличии от неё она является гладкой, дифференцируемой. Потенциал уменьшается, приближаясь к фракталу, и на нём равен 0, и бесконечно увеличивается при бесконечном удалении от фрактала. Потенциал задаёт эквипотенциальные поверхности, которые тем плотнее прилегают к фракталу, чем меньше расстояние до него.

Итак, точка p принадлежит фракталу, если после n итераций она не вышла за пределы сферы радиусом R. Рассмотрим следующую потенциал-функцию:

\begin{aligned}&G(p) = \max\!\left( \log\!\left( \dfrac{\lVert p_n \rVert}{R} \right),\, 0 \right), \\[6pt]&\text{где } G \text{ — монотонная функция, растёт с удалением от фрактала и стремится к 0} \\[4pt]&\text{при приближении к фракталу;} \\[6pt]&\text{при этом } G = 0 \text{ в точках фрактала.}\end{aligned}

Так как G — монотонная функция, то эквипотенциальная поверхность, потенциал точек которой равен нулю (G = 0), приближает поверхность фрактала, причём тем точнее, чем больше итераций мы делаем. Тогда расстояние до фрактала можно аппроксимировать расстоянием до этой поверхности, так как все точки фрактальной поверхности лежат на поверхности уровня 0, а значит и расстояние до неё не больше истинного расстояния до фрактала.

Оценим расстояние до этой поверхности уровня 0.

\begin{aligned}&\text{Пусть } p \text{ — точка на поверхности, ближайшая к нашей точке } p_0 = p + \Delta p. \\[6pt]&\text{Тогда нам нужно оценить величину } d = \lVert \Delta p \rVert, \text{ это и есть наше расстояние.} \\[8pt]&\text{Запишем линейное приближение потенциала в окрестности } p_0: \\[4pt]&G(p) \approx G(p_0) +\left( \dfrac{\partial G}{\partial x} \right)_{p_0} \!\!\Delta x +\left( \dfrac{\partial G}{\partial y} \right)_{p_0} \!\!\Delta y +\left( \dfrac{\partial G}{\partial z} \right)_{p_0} \!\!\Delta z, \\[8pt]&\text{или в более компактном виде:} \\[4pt]&G(p) \approx G(p_0) + \nabla G(p_0) \cdot \Delta p. \\[8pt]&\text{В нашем случае } G(p) = 0, \text{ так как точка } p \text{ лежит на поверхности уровня } 0. \\[4pt]&0 \approx G(p_0) + \nabla G(p_0) \cdot \Delta p. \\[8pt]&\text{Откуда получаем:} \\[4pt]&|G(p_0)| \approx \bigl| \nabla G(p_0) \cdot \Delta p \bigr| = \lVert \nabla G(p_0) \rVert \, \lVert \Delta p \rVert \, |\cos\phi| \le \lVert \nabla G(p_0) \rVert \, \lVert \Delta p \rVert, \\[8pt]&\Rightarrow \lVert \Delta p \rVert \ge \dfrac{|G(p_0)|}{\lVert \nabla G(p_0) \rVert}. \\[10pt]&\text{Значит, расстояние до поверхности уровня 0, а тем более до фрактала, можно оценить снизу как } \\[10pt]& d = \lVert \Delta p \rVert \ge \dfrac{|G(p_0)|}{\lVert \nabla G(p_0) \rVert}\end{aligned}

Мы получили общую формулу оценки расстояния снизу для многих фракталов и не только их.

В нашем случае G(p) = \max\!\left( \log\!\left( \dfrac{ \lVert p_n \rVert }{ R } \right),\, 0 \right), то есть получим такую оценку:

\left| \dfrac{ G }{ \lVert \nabla G \rVert } \right|= \dfrac{ \log\!\left( \dfrac{ \lVert p_n \rVert }{ R } \right) } { \left\lVert \nabla \log\!\left( \dfrac{ \lVert p_n \rVert }{ R } \right) \right\rVert }=\dfrac{ \log\!\left( \dfrac{ \lVert p_n \rVert }{ R } \right) } { \left\lVert \nabla \log\!\left( \lVert p_n \rVert \right) - \nabla \log(R) \right\rVert }=\dfrac{ \log\!\left( \dfrac{ \lVert p_n \rVert }{ R } \right) } { \left\lVert \nabla \log\!\left( \lVert p_n \rVert \right) \right\rVert }.

Обозначим \lVert p_n \rVert = r, тогда:

\dfrac{ \log\!\left( \dfrac{r}{R} \right) } { \lVert \nabla \log r \rVert }=\dfrac{ \log r }{ \lVert \nabla \log r \rVert }-\dfrac{ \log R }{ \lVert \nabla \log r \rVert }.

С ростом n растёт и \lVert \nabla \log r \rVert, поэтому второй член стремится к 0, а остаётся вот этот:

\dfrac{ \log r }{ \lVert \nabla \log r \rVert }=\dfrac{ \log r }{ \dfrac{\partial r}{r} }=\dfrac{ \log r \cdot r }{ \partial r }

Градиент \nabla \log r можно вычислить приближённо или автоматически дифференцировать, например, с помощью дуальных чисел, или просто накопить значение:

\begin{aligned}\partial r_0 &= 1, \\\partial r_1 &= \left( k \, r_0^{k-1} \right) \partial r_0 + 1,\\ \vdots\\\partial r_{i+1} &= \left( k \, r_i^{k-1} \right) \partial r_i + 1.\end{aligned}

Оценка\dfrac{ \log r \cdot r }{ \partial r } хороша вблизи поверхности. При удалении за пределы сферы радиусом R

итерации имеет смысл обрывать и уменьшить оценку, чтобы не перескочить поверхность фрактала. Например, так:

\alpha \cdot \dfrac{ \log r \cdot r }{ \partial r },\quad \text{где } \alpha \log(r) < 1,\quad \alpha = 0.25.

В итоге получаем оценку расстояния:

\text{dist} = 0.25 \cdot \dfrac{ \log r \cdot r }{ \partial r }.

Напишем SDF.

\begin{aligned} & MbulbDeg(\text{vec3}\ p,\ \text{float}\ deg)\ \{ \\[-0.3em] & \quad \text{// получаем сферические координаты} \\ & \quad \text{vec3}\ s = cartezianToSphere(p); \\ & \quad \text{// s = (r, θ, φ)} \\ & \quad \text{float}\ t = \theta \cdot deg \cdot paramX; \\ & \quad \text{float}\ f = \phi \cdot deg \cdot paramY; \\ & \quad \text{// здесь добавлено умножение на параметры paramX, paramY} \\ & \quad \text{// по умолчанию они равны 1, меняя их получим различные формы} \\ & \quad \text{vec3}\ m = r^{deg} \cdot (\, \sin(t)\cos(f),\ \sin(t)\sin(f),\ \cos(t)\,); \\ & \quad \text{return}\ m; \\ & \} \\[1em] & sdMandelBulb(\text{vec3}\ p,\ \text{float}\ deg)\ \{ \\[-0.3em] & \quad \text{vec3}\ p_0 = p; \\ & \quad \text{vec3}\ np = p; \\ & \quad \text{float}\ dr = 1.; \\ & \quad \text{float}\ r; \\ & \quad \text{// цикл по глубине итераций} \\ & \quad \text{for}(\text{int}\ i = 1;\ i < depth;\ i++)\ \{ \\ & \quad\quad r = \lVert np \rVert; \\ & \quad\quad \text{// если расстояние велико — выходим} \\ & \quad\quad \text{if}(r > 5.)\ \text{break}; \\ & \quad\quad np = MbulbDeg(np,\ deg) + p_0; \\ & \quad\quad dr = dr \cdot deg \cdot r^{deg - 1.} + 1.; \\ & \quad \} \\ & \quad \text{float}\ dist = 0.25 \cdot \dfrac{\log(r) \cdot r}{|dr|}; \\ & \quad \text{return}\ dist; \\ & \} \end{aligned}

Получим следующее изображение для степени deg = 8.

Оболочка Мандельброта (степень 8)
Оболочка Мандельброта (степень 8)

Добавим ловушку для раскрашивания, и меняя параметры, получим, например, такие изображения и срезы фракталов.

Аналогично можно построить множества Жюлиа (Juliabulb). Оценку расстояния можно оставить такой же, или сделать небольшое уточнение в производной:

\partial r_0 = 1, \quad\partial r_1 = \left( k \cdot r_0^{\,k-1} \right) \cdot \partial r_0, \quad\dots, \quad\partial r_{i+1} = \left( k \cdot r_i^{\,k-1} \right) \cdot \partial r_i (не добавлять +1, в итерациях).

Получим такие изображения:

Фракталы Мандельброта и Жюлиа на кватернионах

На кватернионах мы можем построить четырёхмерные аналоги множеств Мандельброта и Жюлиа. Визуализировать мы будем, конечно, пересечение фрактала с некоторым трёхмерным подпространством (например, взяв одну из координат равной 0, (x, y, z, 0)).

Для оценки расстояния до фрактала используем формулу, которую мы получили при оценке расстояния до оболочки Мандельброта. Функцию потенциала можно взять ту же. Оценка получится точно такой же:

\text{dist} = 0.25 \cdot \dfrac{ \log r \cdot r }{ \partial r }.

Чтобы написать SDF для четырёхмерного фрактала Мандельброта, нам нужна будет функция возведения в степень кватерниона. Для этого приведём кватернион к тригонометрической форме.

\begin{aligned}q &= a + b i + c j + d k, \quad v = b i + c j + d k, \\[0.3em]r &= \lVert q \rVert, \quad \theta = \arccos\!\left(\dfrac{a}{r}\right), \\[0.3em]q &= r \cdot \left( \cos(\theta),\ \sin(\theta) \cdot \dfrac{v}{\lVert v \rVert} \right), \\[0.3em]q^{n} &= r^{n} \cdot \left( \cos(n \cdot \theta),\ \sin(n \cdot \theta) \cdot \dfrac{v}{\lVert v \rVert} \right)\end{aligned}

Тогда SDF:

\begin{aligned} & qDeg(\text{vec4}\ q,\ \text{float}\ deg)\ \{ \\[-0.3em] & \quad \text{// q = a + b i + c j + d k, v = b i + c j + d k} \\ & \quad \text{float}\ r = \lVert q \rVert; \\ & \quad \text{float}\ \theta = \arccos\!\left(\dfrac{a}{r}\right); \\ & \quad \text{vec4}\ nq = r^{deg} \cdot \left(\cos(\theta),\ \sin(\theta) \cdot \dfrac{v}{\lVert v \rVert}\right); \\ & \quad \text{return}\ nq; \\ & \} \\[1em] & sdMandelQ(\text{vec3}\ p,\ \text{float}\ deg)\ \{ \\[-0.3em] & \quad \text{vec4}\ q = \text{vec4}(p,\ 0); \\ & \quad \text{vec4}\ nq = q; \\ & \quad \text{float}\ dr = 1.;\ \text{float}\ r; \\ & \quad \quad \text{for}(\text{int}\ i = 1;\ i < depth;\ i++)\ \{ \\ & \quad\quad r = \lVert nq \rVert; \\ & \quad\quad \quad\quad \text{if}(r > 2.)\ \text{break}; \\ & \quad\quad nq = qDeg(nq,\ deg) + q; \\ & \quad\quad dr = dr \cdot deg \cdot r^{deg - 1.} + 1.; \\ & \quad \} \\ & \quad \text{float}\ dist = 0.25 \cdot \dfrac{\log(r) \cdot r}{|dr|}; \\ & \quad \text{return}\ dist; \\ & \} \end{aligned}

Получим такие изображения:

Аналогично построим множества Жюлиа на кватернионах.

Гибриды

Рассмотрим теперь естественные способы получения новых фрактальных множеств и оценок расстояний до них, исходя из уже имеющихся. Такие фрактальные множества называют гибридами. Я приведу два способа создания гибридов. В действительности их намного больше, но если их все описывать, то статья вырастет раза в два, поэтому я решил ограничиться двумя, наиболее часто используемыми.

Для примеров нам нужна будет губка Менгера, поэтому сначала найдём её SDF как аттрактора соответствующей IFS (до этого мы уже находили SDF губки исходя из других соображений).

\begin{aligned}& sdMenger(p)\ \{ \\[-0.3em]& \quad \text{float}\ scale = escapeForce + 1.;\ \text{// коэффициент гомотетии} \\& \quad \text{vec3}\ v = \text{vec3}(1,\ 1,\ 1); \\& \quad \text{vec3}\ pt = p; \\& \quad \text{for}(\text{int}\ i = 0;\ i < dep;\ ++i)\ \{ \\& \quad\quad pt = transformation(pt);\ \text{// поворот и перенос точки} \\[0.3em]& \quad\quad \text{if}(pt.x - pt.y < 0.0)\ \{\, pt.xy = pt.yx;\,\}\ \text{// симметрии} \\& \quad\quad \text{if}(pt.x - pt.z < 0.0)\ \{\, pt.xz = pt.zx;\,\} \\& \quad\quad \text{if}(pt.y - pt.z < 0.0)\ \{\, pt.yz = pt.zy;\,\} \\[0.3em]& \quad\quad pt = transformation(pt); \\[0.3em]& \quad\quad pt.xy = scale \cdot pt.xy + (1. - scale) \cdot v.xy; \\& \quad\quad pt.z = scale \cdot pt.z; \\[0.3em]& \quad\quad \text{if}(pt.z > (scale - 1.0) \cdot 0.5 \cdot v.z)\ \{ \\& \quad\quad\quad pt.z = pt.z - (scale - 1.0) \cdot v.z; \\& \quad\quad \} \\& \quad \} \\[0.3em]& \quad \text{return}\ |\lVert pt \rVert - 1.| \cdot scale^{-dep}; \\& \}\end{aligned}

Получим, например, такие изображения губки:

Первый способ: алгебраическая гибридизация.

Пусть есть два фрактала, заданные правилами F и G:

p_{n+1} = F(p_n), \quad \text{и} \quad p_{n+1} = G(p_n)

Тогда, используя алгебраические операции, определённые на точках, построим новый фрактал H. Причём если мы умеем вычислять скалярные производные, градиенты или якобианы преобразований, то и производную алгебраического гибрида легко найдём, используя свойства производной. Количество и тип используемых операций ограничены лишь вашей фантазией. Например, можно взять сумму F и G.

H(p) = \alpha \cdot F(p) + (1 - \alpha) \cdot G(p), \quad \text{где обычно } \alpha \in [0, 1],\\[0.3em]p_{n+1} = H(p_n)

Естественно, можно использовать более двух фрактальных преобразований. Сумма хорошо подходит для морфинга фракталов.

\begin{aligned}& a \in [0, 5.5] \\[0.3em]& t_i = \text{clamp}\Big( 2 \cdot (1 - |a - 1.5 \cdot i|),\ 0,\ 1 \Big), \quad i \in \{0, 1, 2, 3, 4\} \\[0.3em]& \text{// Тут в цикле выполняем вычисления} \\[0.3em]& p = t_0 \cdot \text{tetraSerp}(p) + t_1 \cdot \text{Menger}(p) \\ & \quad + t_2 \cdot \text{oktaedrSerp}(p) + t_3 \cdot \text{icosaedrSerp}(p) + t_4 \cdot \text{dodecaedrSerp}(p); \\[0.3em]& \text{// домножаем производные на коэффициенты и накапливаем в общей производной dr} \\[0.3em]& d_i = d_i \cdot t_i, \quad i \in \{0, 1, 2, 3, 4\} \\[0.3em]& dr = dr \cdot (d_0 + d_1 + d_2 + d_3 + d_4); \\[0.3em]& r = \lVert p \rVert; \\& \text{if } \lVert p \rVert > 5. \text{ break}; \\[0.3em]& \text{// И выходя из цикла, вычисляем расстояние} \\[0.3em]& \text{return } \dfrac{r}{dr};\end{aligned}

Для примера, я записал такое видео: Rutube или Youtube.

Второй способ: композиция преобразований.

Идея в том, чтобы в некоторой последовательности выполнять различные фрактальные преобразования над точкой, а производную вычислять по правилу цепочки.

Например, последовательно выполняя boxFold, инверсию и добавив параллельный перенос и поворот над точками, можно получить так называемые псевдоклейновские фракталы, которые подобны предельным множествам групп Клейна.

А вот, например, какие изображения можно получить, если итерация будет состоять из других комбинаций преобразований.

Три октаэдра Серпинского и две губки Менгера
Три октаэдра Серпинского и две губки Менгера
Губка Менгера, boxFold, sphereFold
Губка Менгера, boxFold, sphereFold

Эти два способа можно также комбинировать.

Заключение

Вот и подошло к завершению наше путешествие в мир трехмерных фрактальных множеств. Очень многое осталось за пределами нашего внимания: фракталы Клейна-Мёбиуса, фракталы Ньютона, странные аттракторы и многое другое. Эти и многие другие визуализации фракталов, ландшафтов, узлов, автоматов и других математических объектов вы можете найти на моём сайте

Надеюсь, вам было интересно. И если так, то ещё увидимся в других статьях :)

Благодарю вас, и удачи! :)

Комментарии (7)


  1. Tzimie
    23.10.2025 14:57
    #29004236

    А что бы такого построить на простых числах или нулях функции Римана?


    1. Andrey-82 Автор
      23.10.2025 14:57
      #29004314

      Вопрос интересный, но до конца не понятный! :)
      С простыми числами можно много чего интересного делать, мне почему-то сразу вспомнилась скатерть Улама. А вот про нули дзета-функции не знаю, можно пробовать саму функцию визуализировать


      1. Tzimie
        23.10.2025 14:57
        #29004626

        Я много пытался, но фракталов не получалось


  1. Real3L0
    23.10.2025 14:57
    #29004604

    Эти фигуры быстро отрисовываются?

    Я к тому, что можно на этом screensaver сделать?


    1. Andrey-82 Автор
      23.10.2025 14:57
      #29004774

      Да, отрисовываются быстро, но зависит от видеокарты. Если видеокарта не хуже, чем моя на ноуте (NVIDIA GeForce RTX 3050 Laptop GPU), то все хорошо будет и можно делать screensaver.


  1. vater_nicht
    23.10.2025 14:57
    #29006872

    На 3-d печать и в серию


    1. Andrey-82 Автор
      23.10.2025 14:57
      #29007310

      :)