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

Меня всё также зовут Андрей Гринблат, и в первой части я начал рассказывать о такой технологии, как ray marching, и о нормированных пространствах. В этой части начнём с построения простых геометрических фракталов — губки Менгера и тетраэдра Серпинского, затем построим IFS-фракталы, рассмотрим технику орбитальных ловушек, и в завершение построим фрактал «Ящик Мандельброта», или Мандельбокс.

Губка Менгера

Реализовывать губку будем в пространстве с предельной нормой.

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

Множество, удаляемое из куба при первой итерации
Множество, удаляемое из куба при первой итерации

SDF операции вычитания объектов было описано ранее, а сам «крест» есть объединение трёх «вытянутых» кубов. Полученное ��наивное» SDF не является точным, но для наших целей более чем удовлетворительно.

SDF креста получаем так: вытягиваем по осям, причём вытянуть нужно более чем в три раза, поэтому знаменатель сделаем равным, например, 4:

p=\left( {x}, y, z \right),\\p_1 = \left( \frac{x}{4}, y, z \right),\;p_2 = \left( x, \frac{y}{4}, z \right),\;p_3 = \left( x, y, \frac{z}{4} \right)\\d_i = sdSphere(p_i, \frac{1}{3}),\; i \in \{1, 2, 3\}sdKrest(p) = \min(d_1, d_2, d_3)

Получим SDF в первой итерации:

sdMenger(p) = \max\bigl( sdSphere(p, 1),\; -sdKrest(p) \bigr)

Отрисуем первую итерацию

Губка Менгера, шаг 1
Губка Менгера, шаг 1

Видим, что внутри губки не учитываются тени и дополнительное рассеивание света: выглядит так, будто освещенность внутри губки, такая же, как и снаружи. Добавим простую ложную окклюзию, это сделает изображение более реалистичным. Для этого умножим вектор цвета фрагмента (RGB-ую часть в данном случае) на что-нибудь, например, такое:

\frac{AO}{\log(i + 100)}, где AO — параметр окклюзии, а i — количество шагов по лучу.

Губка Менгера, шаг 1, + ложная окклюзия
Губка Менгера, шаг 1, + ложная окклюзия

Следующий и остальные шаги в построении губки — это уменьшение «креста» и удаление его из уменьшенных кубов, то есть масштабирование и повторение домена. Будем как бы слоями удалять «кресты»: на первом шаге — слой из 1 «креста», на втором — из 20, уменьшенных втрое, на третьем — из 400, уменьшенных ещё втрое, и так далее Для этого сделаем, например, такой цикл:

\begin{aligned}&\text{float } t; \\[4pt]&\text{float } s = 1.; \\[4pt]&\text{float } d = sdSphere(p, 1.); \\[6pt]&\text{for (int i = 0; i < depth; i++) \{} \\[4pt]&\qquad //\; тут\; depth\;-\;глубина,\; количество\; итераций \\[4pt]&\qquad \text vec3\; ps = p \cdot s; \\[4pt]&\qquad //\; масштабируем,\; раздувае��\; область \\[4pt]&\qquad t = \dfrac{sdKrest\!\left(ps - \left[\dfrac{ps + 1.}{2.}\right] \cdot 2.\right)}{s}; \\[4pt]&\qquad //\; повторяем\; домен\; и\; обратно\; сжимаем\; область \\[4pt]&\qquad s = s \cdot 3.; \\[4pt]&\qquad d = \max(d, -t); \\[4pt]&\text{\}} \\[4pt]&\text{return } d;\end{aligned}

Получили последовательные приближения губки.

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

Простейшие преобразования и операции с фракталом

Описанные ниже преобразования и операции можно применять не только к губке Менгера, но и ко многим другим фракталам. Это позволяет увеличить вариативность получаемых форм. Существует и множество других интересных преобразований: искажения, добавление шума в форму — всё ограничено лишь вашей фантазией!

Но для начала сделаем немного более привлекательным окружающее губку пространство: на плоскости сделаем паркет и добавим текстуру на небесную сферу.

Сделаем срез губки.

Срез губки Менгера плоскостью x + y + z = 0
Срез губки Менгера плоскостью x + y + z = 0

Посмотрим на губку в различным пространствах.

Давайте теперь в SDF губки вместо координат текущей точки будем передавать координаты точки, которая получена из текущей с помощью переноса на некоторый вектор и поворота вокруг выбранной оси на некоторый угол.

С переносами всё очевидно: просто прибавляем вектор переноса. А вот повороты реализуем через умножение кватернионов, которое нам также понадобится для построения трёхмерной проекций, четырёхмерных фракталов Жулиа и Мандельброта. Кроме того, группа кватернионов SU(2), являющаяся двулистным накрытием группы вращений в трёхмерном пространстве SO(3), более удобна в использовании. По сравнению с матрицами из SO(3) кватернионы обеспечивают более точные и стабильные вращения, не накапливая численных ошибок и всегда выполняя поворот по кратчайшему пути.

Такие преобразования с точками приводят к изменениям формы фрактала. Примеры: один, два.

Рассмотрим построение других фракталов. Начнём с тетраэдра Серпинского.

Тетраэдр Серпинского

Нулевая итерация тетраэдра Серпинского — это обычный тетраэдр. Давайте найдём его SDF.

Расстояние от точки до тетраэдра — это минимальное из расстояний от точки до его граней. Если точка находится внутри тетраэдра, то минимальное из расстояний до граней окажется равным минимальному из расстояний до плоскостей, содержащих эти грани. Находим его и берём со знаком «-», так как область внутренняя. Если же точка находится снаружи, то придётся находить ра��стояние именно до граней.

Выясним расположение нашей точки и самого тетраэдра относительно его граней, посмотрим, по одну или по разные стороны они лежат. Для этого вычислим пять определителей размером 4×4.

v_0, v_1, v_2, v_3 — вершины тетраэдра;

P — текущая точка;

V = \begin{vmatrix}v_0 & 1 \\v_1 & 1 \\v_2 & 1 \\v_3 & 1\end{vmatrix}P_0 =\begin{vmatrix}P & 1 \\v_1 & 1 \\v_2 & 1 \\v_3 & 1\end{vmatrix},\;P_1 =\begin{vmatrix}v_0 & 1 \\P & 1 \\v_2 & 1 \\v_3 & 1\end{vmatrix},\;P_2 =\begin{vmatrix}v_0 & 1 \\v_1 & 1 \\P & 1 \\v_3 & 1\end{vmatrix},\;P_3 =\begin{vmatrix}v_0 & 1 \\v_1 & 1 \\v_2 & 1 \\P & 1\end{vmatrix}s_i = \operatorname{sign}(V \cdot P_i)
  • если s_i = -1, то точка и тетраэдр находятся по разные стороны от i-й грани;

  • если s_i = 1, то точка и тетраэдр находятся по одну сторону от i-й грани.

Если все s_i = 1, то текущая точка расположена внутри тетраэдра

Если ровно одно s_i = -1, то в качестве нижней грани расстояния можно взять расстояние от P до плоскости, определяемой вершинами с номерами:

\{0, 1, 2, 3\} \setminus \{i\}

Если ровно два значения отрицательны, i, j, то в качестве нижней грани расстояния можно взять расстояние от P до прямой, определяемой вершинами с номерами:

\{0, 1, 2, 3\} \setminus \{i, j\}

Наконец, если только одно s_i = +1, то в качестве расстояния можно взять расстояние от P до i-ой вершины тетраэдра.

Эти вычисления можно оптимизировать. Например, заменив вычисления определителей на соответствующие вычисления векторных и скалярных произведений, которые также используются в дальнейшем.

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

Нарисуем правильный тетраэдр с вершинам:

V_0 (1, 1, -1), V_1 (1, -1, 1), V_2 (-1, 1, 1), V_3 (-1, -1, -1)
Тетраэдр Серпинского, шаг 0
Тетраэдр Серпинского, шаг 0

Рассмотрим следующую — первую — итерацию тетраэдра Серпинского, она представляет собой тетраэдр, из которого вынули октаэдр. Октаэдр построим как сферу в пространстве с манхэттенской нормой.

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

Тетраэдр Серпинского, шаг 1
Тетраэдр Серпинского, ��аг 1

Немного раскрасим изображение: генерируем текстуры на плоскость и на поверхность тетраэдра, а на небесную сферу натянем фото. И дальше, как вы могли уже догадаться, с помощью операций «повторение домена», масштабирования и удаления будем так же, как и в случае губки Менгера, слоями удалять октаэдры.

\begin{aligned}&\text{float } t; \\[4pt]&\text{float } s = 1.; \\[4pt]&\text{float } d = \text{sdTetra}(p, v_0, v_1, v_2, v_3); \\[6pt]&\text{for (int i = 0; i < depth; i++) \{} \\[4pt]&\qquad //\; тут\; depth\;-\;глубина,\; количество\; итераций \\[4pt]&\qquad \text{vec3 } ps = (p + 1.) \cdot s; \\[4pt]&\qquad //\; масштабируем,\; раздуваем\; область \\[4pt]&\qquad t = \dfrac{\text{sdSphere}\Bigl(ps - \Bigl[\dfrac{ps + 2}{2}\Bigr] \cdot 2 + 1\Bigr)}{s}; \\[4pt]&\qquad //\; повторяем\; домен\; и\; обратно\; сжимаем\; область \\[4pt]&\qquad s = s \cdot 2.; \\[4pt]&\qquad d = \max(d, -t); \\[4pt]&\text{\}} \\[4pt]&\text{return } d;\end{aligned}

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

Для примера рассмотрим ещё фрактал — «многоэтажные дома Дзюбы».

Многоэтажные дома Дзюбы

Для его построения берут куб, из которого удаляют «сферы» (в различных вариантах нормы пространства).

\begin{aligned}&\text{float } t; \\[4pt]&\text{float } s = 1.; \\[4pt]&\text{float } d = \text{sdSphere}(p, 0.67); \\[4pt]&\texttt{//}\;\text{в пространстве с предельной нормой} \\[6pt]&\text{for (int i = 0; i < depth; i++) \{} \\[4pt]&\qquad \texttt{//}\;\text{тут } \texttt{depth}\; \text{- глубина, количество итераций} \\[4pt]&\qquad \text{vec3 } ps = \text{mod}((p \cdot s), 1.) - 0.5; \\[4pt]&\qquad t = \dfrac{\text{sdSphere}\Bigl(p, 0.9 + st\Bigr)}{s}; \\[4pt]&\qquad \texttt{//}\;\text{тут } \texttt{st}\; \text{= } x_0 + y_0 + z_0,\; \text{где } (x_0, y_0, z_0)\; \text{- вектор переноса} \\[4pt]&\qquad s = s \cdot (\text{escapeForce} + 1.); \\[4pt]&\qquad \texttt{//}\;\texttt{escapeForce}\; \text{- параметр масштабирования} \\[4pt]&\qquad d = \max(d, -t); \\[4pt]&\text{\}} \\[4pt]&\text{return } d;\end{aligned}

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

Перейдём к построению октаэдра Серпинского.

Октаэдр Серпинского

Рассмотрим другие, возможно, более естественные методы построения трёхмерных фракталов.

Октаэдр Серпинского представляет собой притягивающее множество (аттрактор) системы итерируемых функций (IFS). Эта система, в данном случае, состоит из сжимающих гомотетий с центрами в вершинах октаэдра и коэффициентом сжатия, равным 0,5.

Будем искать расстояние от точки P до аттрактора последовательными приближениями. Октаэдр Серпинского состоит из шести, уменьшенных в два раза, его копий. Выберем из них ближайшую к нашей точке. Для этого найдём расстояния до вершин исходного октаэдра и возьмём ту копию, которой принадлежит вершина с наименьшим расстоянием до неё. Задача свелась к нахождению расстояния от точки P до этой копии октаэдра. Повторяя процесс, в пределе получим последовательность расстояний, сходящуюся к искомому расстоянию. Кроме того, каждую такую последовательность можно представить как последовательность чисел из {0, 1, 2, 3, 4, 5} и использовать это, например для определения цвета фрагмента. Напишем почти точную (точную в пределе) SDF октаэдра Серпинского.

\begin{aligned}&\text{float } \text{sdOktSerp}(p)\; \{ \\[4pt]&\qquad \text{float } s = \text{sdSphere}_{manhattan}(p, 1.1); \\[4pt]&\qquad \texttt{// оптимизируем вычисления} \\[4pt]&\qquad \text{if } (s > 0.)\; \text{return } s + 0.05; \\[4pt]&\qquad \texttt{// первый шаг делаем сразу к октаэдру} \\[6pt]&\qquad \text{float } \text{scale} = 1. / \text{escapeForce}; \\[4pt]&\qquad \texttt{// коэффициент гомотетии} \\[4pt]&\qquad \text{vec3 } v[6]; \\[4pt]&\qquad \texttt{// массив вершин октаэдра} \\[4pt]&\qquad v_{0} = \text{vec3}(1, 0, 0); \\[4pt]&\qquad v_{1} = \text{vec3}(0, 1, 0); \\[4pt]&\qquad v_{2} = \text{vec3}(0, 0, 1); \\[4pt]&\qquad v_{3} = -v_{0}; \\[4pt]&\qquad v_{4} = -v_{1}; \\[4pt]&\qquad v_{5} = -v_{2}; \\[6pt]&\qquad \text{vec3 } \text{curV}; \\[4pt]&\qquad \text{float } \text{dist}[6]; \\[4pt]&\qquad \text{float } d = 100.; \\[6pt]&\qquad \text{for (int i = 0; i < 15; ++i) \{} \\[4pt]&\qquad\qquad \text{for (int j = 0; j < 6; ++j) \{} \\[4pt]&\qquad\qquad\qquad \texttt{// находим ближайшую вершину} \\[4pt]&\qquad\qquad\qquad \text{dist}_{j} = \lVert p - v_{j} \rVert; \\[4pt]&\qquad\qquad\qquad \text{if } (d > \text{dist}_{j})\; \{ \\[4pt]&\qquad\qquad\qquad\qquad d = \text{dist}_{j}; \\[4pt]&\qquad\qquad\qquad\qquad \text{curV} = v_{j}; \\[4pt]&\qquad\qquad\qquad \} \\[4pt]&\qquad\qquad \} \\[6pt]&\qquad\qquad v_{0} = \text{scale} \cdot v_{0} + (1. - \text{scale}) \cdot \text{curV}; \\[4pt]&\qquad\qquad v_{1} = \text{scale} \cdot v_{1} + (1. - \text{scale}) \cdot \text{curV}; \\[4pt]&\qquad\qquad v_{2} = \text{scale} \cdot v_{2} + (1. - \text{scale}) \cdot \text{curV}; \\[4pt]&\qquad\qquad v_{3} = \text{scale} \cdot v_{3} + (1. - \text{scale}) \cdot \text{curV}; \\[4pt]&\qquad\qquad v_{4} = \text{scale} \cdot v_{4} + (1. - \text{scale}) \cdot \text{curV}; \\[4pt]&\qquad\qquad v_{5} = \text{scale} \cdot v_{5} + (1. - \text{scale}) \cdot \text{curV}; \\[4pt]&\qquad \} \\[6pt]&\qquad \text{return } d; \\[4pt]&\}\end{aligned}

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

Октаэдр Серпинского
Октаэдр Серпинского

Рассмотрим теперь другой, более эффективный метод построения IFS трёхмерных фракталов.

Длина пробега точки

Для простоты пока ограничимся рассмотрением IFS, каждая функция которой обратима. В случае октаэдра имеем систему гомотетий с теми же, что и ранее, центрами, но коэффициентами, равными 2. Для такой IFS октаэдр Серпинского будет отталкивающим множеством (репеллером).

Если точка принадлежит октаэдру, то её образ при действии гомотетии с центром ближайшим к ней также будет принадлежать октаэдру. Всякая другая точка будет удаляться от репеллера.

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

На рисунке изображена, в некотором приближении, салфетка Серпинского и точка P, расстояние от которой до фрактала мы и хотим узнать. Ближайшая вершина для P — это вершина 1, поэтому действуем гомотетией 1. Точка переходит в красную точку 1, для неё ближайшая вершина 0, под действием гомотетии с центром в 0 точка переходит в точку 10; затем ближайшая вершина 2, точка 10 переходит в точку 102, снова ближайшая вершина 2, 102 переходит в 1022, и так далее.

Очевидно, что каждый раз расстояние до фрактала (длины красных отрезков) увеличивается вдвое. Таким образом, расстояние от точки 10 222 до фрактала в 25 раз больше исходного искомого расстояния. Тогда можно оценить снизу это расстояние.

\begin{aligned}&\text{Пусть } d \text{ — искомое расстояние от точки } P \text{ до фрактала,} \\[4pt]&\quad P_{n} \text{ — энная итерация точки } P, \\[4pt]&\quad d_{n} \text{ — расстояние от точки } P_{n} \text{ до фрактала,} \\[4pt]&\quad Z \text{ — произвольная точка фрактала или некоторая точка во «внутренней» области фрактального множества,} \\[4pt]&\quad s \text{ — коэффициент гомотетии } (s > 1), \\[4pt]&\quad \text{diam — диаметр фрактального множества.} \\[10pt]&\text{Тогда:} \\[6pt]&\quad d_{n} = d \cdot s^{n} \;\Rightarrow\;d = \dfrac{d_{n}}{s^{n}} \;\ge\;\dfrac{\lVert Z - P_{n} \rVert - \text{diam}}{s^{n}}, \\[10pt]&\text{и кроме того:} \\[6pt]&\quad d = \lim_{n \to \infty} \dfrac{\lVert Z - P_{n} \rVert - \text{diam}}{s^{n}}.\end{aligned}

В качестве точки Z можно взять любую из вершин октаэдра, а в качестве диаметра — расстояние между противоположными вершинами.

Коэффициент гомотетий s = 2.

В нашем случае в качестве Z можно выбрать начало координат, и в качестве diam взять половину диаметра = 1.

Таким образом, получим:

\begin{aligned}&\mathrm{sdOktSerp}(P) = \lim_{n \to \infty} \dfrac{\lVert P_{n} \rVert - 1}{2^{n}}, \\[10pt]&\text{и для любого конечного натурального}\ n = k, \\[4pt]&\mathrm{ошибка\ будет\ равна}\ \dfrac{1}{2^{k}}.\end{aligned}

Так, например, если мы хотим иметь величину ошибки не больше, чем 0,00001, то нужно проитерировать точку 17 или более раз.

Для оптимизации вычислений и для получения большего разнообразия форм воспользуемся симметриями октаэдра. Идея состоит в том, что вместо того, чтобы вычислять расстояния до вершин октаэдра и находить меньшее из них, иметь дело всегда только с одной выбранной вершиной (например, (1, 0, 0)). Для этого всякий раз, когда ближайшая вершина не (1, 0, 0), симметрично отражать точку P так, чтобы ближайшая к ней вершина отразилась в (1, 0, 0).

\begin{aligned}&\text{float sdOktSerp}(p) \; \{ \\[4pt]&\qquad \text{float } \text{scl} = \text{escapeForce}; \;\texttt{// коэффициент гомотетии} \\[4pt]&\qquad \text{vec3 } v = \text{vec3}(1, 0, 0); \;\texttt{// выбранная вершина октаэдра} \\[4pt]&\qquad \text{vec3 } pt = p; \\[4pt]&\qquad \text{for (int i = 0; i < dep; ++i) \{} \\[4pt]&\qquad\qquad pt = \text{transformation}(pt); \;\texttt{// поворот и перенос точки} \\[4pt]&\qquad\qquad \text{if } (pt.x + pt.y < 0.)\; \{ pt.xy = -pt.yx; \} \;\texttt{// симметрии} \\[4pt]&\qquad\qquad \text{if } (pt.x + pt.z < 0.)\; \{ pt.xz = -pt.zx; \} \\[4pt]&\qquad\qquad \text{if } (pt.x - pt.y < 0.)\; \{ pt.xy = pt.yx; \} \\[4pt]&\qquad\qquad \text{if } (pt.x - pt.z < 0.)\; \{ pt.xz = pt.zx; \} \\[4pt]&\qquad\qquad pt = \text{scale} \cdot pt + (1. - \text{scale}) \cdot v; \\[4pt]&\qquad \} \\[4pt]&\qquad \text{return } \lvert \lVert pt \rVert - 1. \rvert \cdot \text{scale}^{-dep}; \\[4pt]&\}\end{aligned}

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

Икосаэдр Серпинского

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

\begin{aligned}&\phi = \dfrac{\sqrt{5.0} + 1.0}{2.0}; \;\texttt{// число Фибоначчи} \\[6pt]&\text{vec3 } n_{1} = \text{normalize}(-\phi,\ \phi - 1,\ 1); \;\texttt{// нормальные вектора} \\[4pt]&\text{vec3 } n_{2} = \text{normalize}(1,\ -\phi,\ \phi + 1); \;\texttt{// плоскостей симметрии} \\[4pt]&\text{vec3 } n_{3} = \text{normalize}(0,\ 0,\ -1); \\[8pt]&\text{float sdIcosahSerp}(p)\; \{ \\[4pt]&\qquad \text{float } \text{scl} = \text{escapeForce}; \;\texttt{// коэффициент гомотетии} \\[4pt]&\qquad \text{vec3 } v = \text{vec3}(0.8507,\ 0.5257,\ 0); \;\texttt{// выбранная вершина икосаэдра} \\[4pt]&\qquad \text{vec3 } pt = \lvert p \rvert; \;\texttt{// отражаем точку P в положительную область} \\[6pt]&\qquad pt = pt - 2.0 \cdot \max(0.0,\ \text{dot}(pt,\ n_{1})) \cdot n_{1}; \;\texttt{// симметрии} \\[4pt]&\qquad pt = pt - 2.0 \cdot \max(0.0,\ \text{dot}(pt,\ n_{2})) \cdot n_{2}; \\[4pt]&\qquad pt = pt - 2.0 \cdot \max(0.0,\ \text{dot}(pt,\ n_{3})) \cdot n_{3}; \\[4pt]&\qquad pt = pt - 2.0 \cdot \max(0.0,\ \text{dot}(pt,\ n_{2})) \cdot n_{2}; \\[6pt]&\qquad \text{for (int i = 0; i < dep; ++i) \{} \\[4pt]&\qquad\qquad pt = \text{transformation}(pt); \;\texttt{// поворот и перенос точки} \\[4pt]&\qquad\qquad pt = \lvert pt \rvert; \\[4pt]&\qquad\qquad pt = pt - 2.0 \cdot \max(0.0,\ \text{dot}(pt,\ n_{1})) \cdot n_{1}; \\[4pt]&\qquad\qquad pt = \text{scale} \cdot pt + (1. - \text{scale}) \cdot v; \\[4pt]&\qquad \} \\[6pt]&\qquad \text{return } \lvert \lVert pt \rVert - 2. \rvert \cdot \text{scale}^{-dep}; \\[4pt]&\}\end{aligned}

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

Додекаэдр Серпинского

Полностью аналогично напишем SDF для додекаэдра.

\begin{aligned}&\phi = \dfrac{\sqrt{5.0} + 1.0}{2.0}; \;\texttt{// число Фибоначчи} \\[6pt]&\text{vec3 } n_{1} = \text{normalize}(-1,\ \phi - 1,\ \dfrac{1}{\phi - 1}); \;\texttt{// нормальные вектора} \\[4pt]&\text{vec3 } n_{2} = \text{normalize}(\phi - 1,\ \dfrac{1}{\phi - 1},\ -1); \;\texttt{// плоскостей симметрии} \\[4pt]&\text{vec3 } n_{3} = \text{normalize}(\dfrac{1}{\phi - 1},\ -1,\ \phi - 1); \\[8pt]&\text{float sdDodecahSerp}(p)\; \{ \\[4pt]&\qquad \text{float } \text{scl} = \text{escapeForce}; \;\texttt{// коэффициент гомотетии} \\[4pt]&\qquad \text{vec3 } v = \text{vec3}(1,\ 1,\ 1); \;\texttt{// выбранная вершина додекаэдра} \\[4pt]&\qquad \text{vec3 } pt = p; \\[6pt]&\qquad pt = pt - 2.0 \cdot \max(0.0,\ \text{dot}(pt,\ n_{1})) \cdot n_{1}; \;\texttt{// симметрии} \\[4pt]&\qquad pt = pt - 2.0 \cdot \max(0.0,\ \text{dot}(pt,\ n_{2})) \cdot n_{2}; \\[4pt]&\qquad pt = pt - 2.0 \cdot \max(0.0,\ \text{dot}(pt,\ n_{3})) \cdot n_{3}; \\[4pt]&\qquad pt = pt - 2.0 \cdot \max(0.0,\ \text{dot}(pt,\ n_{2})) \cdot n_{2}; \\[6pt]&\qquad \text{for (int i = 0; i < dep; ++i) \{} \\[4pt]&\qquad\qquad pt = \text{transformation}(pt); \;\texttt{// поворот и перенос точки} \\[4pt]&\qquad\qquad \text{for (int j = 0; j < 3; ++j) \{} \\[4pt]&\qquad\qquad\qquad pt = pt - 2.0 \cdot \max(0.0,\ \text{dot}(pt,\ n_{1})) \cdot n_{1}; \\[4pt]&\qquad\qquad\qquad pt = pt - 2.0 \cdot \max(0.0,\ \text{dot}(pt,\ n_{2})) \cdot n_{2}; \\[4pt]&\qquad\qquad\qquad pt = pt - 2.0 \cdot \max(0.0,\ \text{dot}(pt,\ n_{3})) \cdot n_{3}; \\[4pt]&\qquad\qquad \} \\[4pt]&\qquad\qquad pt = \text{scale} \cdot pt + (1. - \text{scale}) \cdot v; \\[4pt]&\qquad \} \\[6pt]&\qquad \text{return } \lvert \lVert pt \rVert - 2. \rvert \cdot \text{scale}^{-dep}; \\[4pt]&\}\end{aligned}

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

Орбитальные ловушки

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

Под орбитой точки x_0 будем понимать множество точек f(x_0), f(f(x_0)), f(f(f(x_0))), …, f^{(n)}(x_0), …, где f — функция, действующая на пространстве. В примерах выше f — это гомотетия в выбранной вершине, а точнее гомотетия вместе с симметриями, предварительно применяемыми к точке. Также. если мы не использовали симметрии, то на каждой итерации f — это гомотетия с центром ближайшим к точке.

Идея орбитальных ловушек состоит в следующем. Выбираем какое‑то множество в качестве ловушки (сферу, плоскость, точку или даже шум, фрактал и так далее) и фиксируем, насколько близко к ловушке проходит орбита точки. Дале, в зависимости от этого расстояния, красим точку. Также мы можем не просто красить точку фрактала, а считать её принадлежащей нашему объекту, даже если она не лежит на фрактале.

Для примера рассмотрим IFS тетраэдра Серпинского:

\begin{aligned}&\text{float sdTetraSerp}(p)\; \{ \\[4pt]&\qquad \text{float } \text{scl} = \text{escapeForce}; \;\texttt{// коэффициент гомотетии} \\[4pt]&\qquad \text{vec3 } v = \text{vec3}(1,\ 1,\ 1); \;\texttt{// выбранная вершина тетраэдра} \\[4pt]&\qquad \text{vec3 } pt = p; \\[6pt]&\qquad \text{for (int i = 0; i < dep; ++i) \{} \\[4pt]&\qquad\qquad pt = \text{transformation}(pt); \;\texttt{// поворот и перенос точки} \\[4pt]&\qquad\qquad \text{if (pt.x + pt.y < 0.) \{ pt.xy = -pt.yx; \}} \;\texttt{// симметрии} \\[4pt]&\qquad\qquad \text{if (pt.x + pt.z < 0.) \{ pt.xz = -pt.zx; \}} \\[4pt]&\qquad\qquad \text{if (pt.y + pt.z < 0.) \{ pt.yz = -pt.zy; \}} \\[4pt]&\qquad\qquad pt = \text{scale} \cdot pt + (1. - \text{scale}) \cdot v; \\[4pt]&\qquad \} \\[6pt]&\qquad \text{return } \lvert \lVert pt \rVert - 1. \rvert \cdot \text{scale}^{-dep}; \\[4pt]&\}\end{aligned}

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

Тетраэдр Серпинского
Тетраэдр Серпинского

Добавим теперь орбитальную ловушку, например, шар с центром в выбранной вершине. Для этого в цикл добавим:

\text{orbital} = \min\!\bigl(\text{orbital},\; \text{sdSphere}(pt - v,\; 0.2)\bigr);

Возвращать будем не просто d — расстояние до фрактала, а расстояние до фрактала с повторяющимися образами ловушки.

\begin{aligned}&d = \lvert \lVert pt \rVert - 1. \rvert \cdot \text{scale}^{-dep}; \\[6pt]&\text{return } \min\!\left(d,\; \dfrac{\text{orbital}}{5}\right);\end{aligned}

SDF с ловушкой не является точным, для большей точности его можно поделить на градиент поля, или уменьшить шаг сэмпла луча, или просто подобрать скаляр, на который поделить расстояние от орбиты до ловушки (в данном примере с тетраэдром поделить на 5).

Получим такое изображение для евклидовой нормы. Видим, как динамика системы гомотетий повторяет образ ловушки на разных уровнях.

Тетраэдр Серпинского с орбитальной ловушкой
Тетраэдр Серпинского с орбитальной ловушкой

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

Давайте сделаем какие-нибудь преобразования над тетраэдром. Например, можем получить такое изображение:

Дерево
Дерево

А теперь то же самое, но с орбитальной ловушкой.

Дерево с плодами
Дерево с плодами

Добавим цвет ловушке и фракталу.

Дерево с плодами (раскрашенное)
Дерево с плодами (раскрашенное)

Мандельбокс

В феврале 2010 года Том Лоу (Tom Lowe) представил найденный им фрактал — коробку Мандельброта. Этот фрактал, как и фрактал Мандельброта, является картой для соответствующих множеств Жулиа, но, в отличие от множества Мандельброта, естественным образом определяется для любого количества измерений.

Точками мандельбокса являются точки, орбиты которых не уходят на бесконечность при итерировании функции:

p_{n+1} = \text{scale} \cdot \text{spherefold}\!\bigl(\text{boxfold}(p_n)\bigr) + p_0 \qquad (*)

где:

  • scale — масштабный фактор, гомотетия с центром в 0 (все точки на рисунке ниже);

  • boxfold — это симметрии относительно сторон коробки, если точка снаружи (серые точки), и тождественное преобразование, если она внутри коробки (не серые);

  • spherefold — это гомотетия с центром в 0 и коэффициентом k^2 для точек, лежащих внутри сферы, длина радиус-вектора которых не больше, чем k (красные точки), и инверсия, если длина радиус-вектора ≥ k и ≤ L (зелёные точки).

Рисунок для простоты и наглядности показывает двумерный случай, в трёхмерном всё аналогично.

L, k, R (радиус большей сферы, на рисунке R = L) и scale — параметры для фрактала.

Давайте теперь найдём SDF мандельбокса.

По построению очевидно, что мандельбокс симметричен относительно ядра (куба, сферы и 0). Для примера построим двумерный мандельбокс.

Пусть точка m — это точка фрактала, ближайшая к данной точке P . Пусть P_n и m_n — образы точек P и m после n‑кратного применения преобразования (*) соответственно. Давайте оценим, как с каждой итерацией будет меняться расстояние между P_i и m_i. Если точка находится за пределами фрактального множества (за областью 3·scale·L), то в каждой итерации расстояние будет увеличиваться примерно в scale раз, в пределе ровно в scale раз.

Если точка внутри области, то при симметриях расстояние не будет меняться, при гомотетии будет увеличиваться в k^2 раз, при инверсии будет увеличиваться примерно в L^2/|Pi|^2 (так как точки достаточно близки), и наконец, добавление первоначальных точек, увеличит расстояние на 1. То есть расстояние увеличится в df_n раз, где df_n — скалярная производная нашего преобразования, которая неограниченно растёт с ростом n.

После применения итераций, так как m — точка фрактала, то и m_n — точка фрактала, тогда получим, например, такую картину:

\begin{aligned}&\text{Пусть } d \text{ — диаметр фрактала, } \quad x = \lVert P - m \rVert \text{ — искомое расстояние, тогда:} \\[6pt]&\lVert P_n \rVert - \dfrac{d}{2} \le \lVert P_n - m_n \rVert, \\[6pt]&\dfrac{\lVert P_n \rVert}{df_n} - \dfrac{d}{2\,df_n} \le \dfrac{\lVert P_n - m_n \rVert}{df_n} = x, \\[6pt]&\text{и в пределе получим:} \\[6pt]&\lim_{n \to \infty}\left(\dfrac{\lVert P_n \rVert}{df_n} - \dfrac{d}{2\,df_n}\right)=\lim_{n \to \infty} \dfrac{\lVert P_n \rVert}{df_n}-\lim_{n \to \infty} \dfrac{d}{2\,df_n}= x - 0 = x.\end{aligned}

Напишем функцию расстояния:

\begin{aligned}&\text{float } df = 1.; \\[4pt]&\text{vec3 boxFold}(\text{vec3 } p,\ \text{float } L) \; \{ \\[2pt]&\qquad \text{vec3 } np = p; \\[2pt]&\qquad np = 2 \cdot \text{clamp}(p, -L, L) - p; \\[2pt]&\qquad \text{return } np; \\[2pt]&\}\end{aligned}\begin{aligned} &\text{vec3 sphereFold}(\text{vec3 } p, \text{float } k, \text{float } R) \{ \\[2pt] &\qquad \text{float } pp = \lVert p \rVert \cdot \lVert p \rVert; \\[2pt] &\qquad \text{float } d = R \cdot R; \\[2pt] &\qquad \text{float } rr = k \cdot k; \\[2pt] &\qquad \text{float } t = 1.; \\[2pt] &\qquad \text{if } (pp < rr) \{ \\[2pt] &\qquad\qquad t = \dfrac{d}{rr}; \; df = df \cdot t; \\[2pt] &\qquad\qquad \text{return } p \cdot t; \\[2pt] &\qquad \} \text{ else if } (pp < d) \{ \\[2pt] &\qquad\qquad t = \dfrac{d}{pp}; \; df = df \cdot t; \\[2pt] &\qquad\qquad \text{return } p \cdot t; \\[2pt] &\qquad \} \\[2pt] &\qquad \text{return } p; \\[2pt] &\} \\[4pt]  &\text{vec3 composFolds}(\text{vec3 } p) \{ \\[2pt] &\qquad \text{vec3 } np; \\[2pt] &\qquad np = \text{boxFold}(p); \\[2pt] &\qquad np = \text{sphereFold}(np); \\[2pt] &\qquad \text{return } np; \\[2pt] &\} \\[4pt]  &\text{float sdMandelBox}(\text{vec3 } p, \text{float } \text{scale}, \text{float } L, \text{float } k, \text{float } R) \{ \\[2pt] &\qquad df = 1.; \\[2pt] &\qquad \text{vec3 } np = p; \\[2pt] &\qquad \text{for (int i = 0; i < dep; ++i) \{} \\[2pt] &\qquad\qquad np = \text{composFolds}(np) \cdot \text{scale} + p; \\[2pt] &\qquad\qquad df = df \cdot \lvert \text{scale} + 1. \rvert; \\[2pt] &\qquad\qquad \text{if } (\lVert np \rVert > 100.) \text{ break}; \\[2pt] &\qquad \} \\[2pt] &\qquad \text{return } \dfrac{\lVert np \rVert}{df}; \\[2pt] &\} \end{aligned}

Получим изображение для значений параметров scale = 2, L = 1, k = 0,5, R = 1.

Вот что, например, можно увидеть внутри ящика:

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

Аналогично множеству Мандельброта можно построить множества Джулиа — Джулиабокс, оценку расстояния для этих множеств можно оставить такую же.

В следующей, заключительной части расскажу про построение алгебраических фракталов: оболочки Мандельброта, фракталов Джулиа и Мандельброта на кватернионах и о построении гибридных фракталов.

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

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


  1. goldexer
    14.10.2025 13:47

    Следующим этапом интересно увидеть масштабирование, имитирующее пролёт внутрь, где фрактальная составляющая снова масштабируется и так далее


    1. Andrey-82 Автор
      14.10.2025 13:47

      Для так построенных фракталов, мы можем легко пролетать внутрь него :)


      1. VBDUnit
        14.10.2025 13:47

        Я когда‑то гонял у себя софт Mandelbulb 3D, который позволет летать внутри 3D фракталов. Не всё так просто. Точность. Нам надо лететь внутри всё медленнее и всё точнее, и со все большей точностью различать всё меньшие детали — мы же во всё меньшие и меньшие закоулки закоулков закоулков фрактала забираемся.

        И даже если положение камеры задаётся числами float64, мы, забираясь внутрь такого фрактала, очень быстро исчерпаем точность своего местоположения. И то что это числа с плавающей запятой нам не поможет. Камера начнёт дёргаться, из поверхностей полезут зубчатые артефакты. Бесконечно углубляться не получится.

        Я думаю, тут надо юзать дроби на длинной арифметике. А это существенно всё усложняет, если мы хотим добиться вменяемой скорости.


        1. Andrey-82 Автор
          14.10.2025 13:47

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


    1. Galkihtuw
      14.10.2025 13:47

      Для видеоролика это решается техниками перебазирования системы координат. Периодически сбрасываешь координаты камеры в ноль, а мир двигаешь под ней. Но для интерактивного реалтайма это уже сложнее


  1. CyberWish
    14.10.2025 13:47

    Очень мощная работа!


    1. Andrey-82 Автор
      14.10.2025 13:47

      Спасибо :)


  1. xtraroman
    14.10.2025 13:47

    Красота. А можете порекомендовать проект чтобы рендерить 3д фракталы?


    1. Andrey-82 Автор
      14.10.2025 13:47

      Вообще, говоря есть около 5, на мой взгляд, неплохих приложений для рендера 3D фракталов. Плюс минус вот эти: Apophysis 3D, ChaosPro, Mandelbulb 3D, Mandelbulber, Fragmentarium. Если у вас видео карта не хуже, чем моя на ноуте (NVIDIA GeForce RTX 3050 Laptop GPU), то скорее всего всё будет отлично, с любой из этих программ.

      Также, я веду свой проект, он в целом посвящён математике, физике, а не только фракталам, но там можно рендерить фракталы :). Вот ссылка https://spherus.ru/
      Проект ещё в стадии разработки, нужны оптимизации и тд, и тп (например, может долго грузить при первой загрузки 3D сцены).


  1. Galkihtuw
    14.10.2025 13:47

    Ray marching для фракталов это чистое искусство. SDF позволяет описывать безумно сложные формы парой строк кода, но за этой простотой скрывается очень нетривиальная математика

    Отдельное спасибо за разбор методов оценки расстояния, это самая суть


    1. Andrey-82 Автор
      14.10.2025 13:47

      Спасибо вам, за ваш интерес! :)
      Полностью согласен, что это во многом искусство.


  1. Jijiki
    14.10.2025 13:47

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

    мне нравятся эффекты постобработки как в фильме джон уик там вообще и свет и обьемный туман тюнили как я понял


    1. Andrey-82 Автор
      14.10.2025 13:47

      Не совсем понял, как объёмный туман, блум, сглаживание и пр, противоречит raymarching? :)