Фракталы, построенные через итерируемые системы функций (Iterated Function Systems, IFS), были формализованы Майклом Барнсли в 1980-х годах. Классический пример — папоротник Барнсли: несколько простых аффинных преобразований, применяемых случайным образом миллионы раз, порождают сложный самоподобный узор.

В 1992 году Скотт Дрейвс (Scott Draves) опубликовал под лицензией GPL первую полную реализацию fractal flame algorithm — расширения классического IFS, которое он разрабатывал с 1991 года, будучи стажёром в NTT-Data (Токио), опираясь на более ранние эксперименты с двумерными IFS, начатые ещё в 1987 году в Brown University. Алгоритм лёг в основу программы Apophysis (2004, Марк Таунсенд) и — позже, в 1999 году, — самостоятельного проекта-продолжения того же Дрейвса, распределённого скринсейвера Electric Sheep. Он сохраняет ядро классического IFS (случайные итерации точки), но добавляет три существенных новшества: нелинейные вариации поверх аффинных преобразований — источник органических, «пламенных» форм вместо чисто геометрических; накопление плотности и цвета в буфере-гистограмме вместо бинарной отрисовки точек; логарифмическое тональное отображение при рендеринге, которое делает видимыми одновременно и плотные, и разреженные области изображения. Код, представленный в примере к лекции, реализует ровно этот алгоритм.

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

Отображение $f: \mathbb{R}^2 \to \mathbb{R}^2$ называется сжимающим (contraction), если существует константа $0 \le s < 1$ такая, что для любых точек $p, q$:
$$| f(p) - f(q) | \le s \cdot | p - q |$$
По теореме Банаха о неподвижной точке, у каждого сжимающего отображения есть единственная неподвижная точка $p^*$, и последовательность $p, f(p), f(f(p)), \dots$ сходится к ней независимо от выбора начальной точки $p$.
Итерируемая система функций — это конечный набор сжимающих отображений $\{F_1, F_2, \dots, F_n\}$ вместе с вероятностями $\{p_1, \dots, p_n\}$, $\sum p_i = 1$. Ключевой факт теории IFS: у такой системы существует единственное непустое компактное множество $A$ (аттрактор), инвариантное относительно объединения отображений:
$$A = \bigcup_{i=1}^{n} F_i(A)$$
Аттрактор — это и есть тот самый фрактальный узор, который мы хотим получить.
Прямое построение $A$ через объединения множеств вычислительно неудобно. Барнсли предложил стохастический алгоритм, который приближает $A$ облаком точек:

1. Взять произвольную точку $x_0$ в двусторонне-единичном квадрате (bi-unit square), то есть $x_0, y_0 \in [-1, 1]$ — именно так стартует и оригинальный псевдокод, и разбираемая реализация (Math.random()*2-1).
2. На шаге $k$ случайно выбрать индекс $i$ с вероятностью $p_i$ и положить $x_{k+1} = F_i(x_k)$.
3. Отметить точку $x_{k+1}$ на изображении.
4. Повторить миллионы раз.

Из-за сжимаемости отображений точка $x_k$ экспоненциально быстро «забывает» начальное положение $x_0$ и с ростом $k$ оказывается сколь угодно близко к аттрактору $A$. Количественно: при типичном коэффициенте сжатия 0.5 уже после 20 итераций точка находится на расстоянии порядка $10^{-6}$ от аттрактора — куда меньше ширины одного пикселя. Отсюда и требование пропускать первые несколько итераций (в оригинальном псевдокоде — именно 20; в разбираемой реализации, с запасом, — 50).

график первых ~10, ~100, ~10000 точек одной и той же chaos game — видно, как облако точек стягивается к финальной форме аттрактора.

В качестве $F_i$ в fractal flame используются аффинные отображения — они дёшевы в вычислении и при этом достаточно выразительны:
$$F_i(x, y) = \begin{pmatrix}a_i & b_i \\ d_i & e_i \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} + \begin{pmatrix} c_i \\ f_i \end{pmatrix} = \begin{pmatrix} a_i x + b_i y + c_i \\ d_i x + e_i y + f_i \end{pmatrix}$$
Геометрически такое отображение — это композиция линейного преобразования (поворот, масштабирование, сдвиг по осям, отражение — всё зависит от матрицы $2\times2$) и параллельного переноса $(c_i, f_i)$.
Условие сходимости. Чтобы хаотическая игра не разбегалась на бесконечность, преобразования должны быть сжимающими в среднем (не обязательно каждое по отдельности). Грубая практическая оценка контрактивности линейной части — сумма модулей коэффициентов $|a|+|b|+|d|+|e|$: если она заметно больше 2, отображение может растягивать плоскость. Строго доказать сжимаемость системы с нелинейными вариациями (о них — ниже) практически невозможно; авторы алгоритма прямо говорят: остаётся полагаться на численно устойчивую реализацию и просто принять, что одни случайные наборы параметров дают удачные узоры, а другие — вырожденные, и предсказать это заранее не получится. В реализации нормировка учтена явно — коэффициенты генерируются случайно, а затем нормируются так, чтобы эта сумма не превышала фиксированный порог:
$$\text{scale} = \frac{|a| + |b| + |d| + |e|}{2}, \qquad \text{если } \text{scale} > 0.85: \quad a, b, d, e \mathrel{{*}{=}} \frac{0.85}{\text{scale}}$$
Коэффициенты сдвига $c, f$ на сходимость не влияют и подбираются свободно (в реализации — в диапазоне, дополнительно уменьшенном множителем 0.6, чтобы преобразования не «разлетались» за пределы холста).
Если ограничиться только аффинными $F_i$, получаются классические IFS-фракталы (папоротник, треугольник Серпинского) — красивые, но геометрически «жёсткие». Идея Дрейвса: после аффинного шага применить дополнительно нелинейную функцию вариации $V$:
$$F_i(x, y) = V_j\big(a_i x + b_i y + c_i,\ d_i x + e_i y + f_i\big)$$
(индекс $j$ у вариации намеренно отделён от индекса $i$ у функции — в общем случае разные функции системы могут использовать одну и ту же или разные вариации). Вариация — произвольная (не обязательно сжимающая) функция $V: \mathbb{R}^2 \to \mathbb{R}^2$, часто определяемая через полярные координаты точки $(x, y)$:
$$r = \sqrt{x^2+y^2}, \qquad \theta = \operatorname{atan2}(x, y)$$
Оговорка по обозначениям: в приложении статьи Дрейвса эта величина формально записана как $\arctan(x/y)$, что само по себе разрывно на границах четвертей; во всех реальных реализациях (включая разбираемую) она вычисляется как двухаргументный $\operatorname{atan2}(x, y)$ — обратите внимание на нестандартный, «перевёрнутый» по сравнению с привычным $\operatorname{atan2}(y, x)$ порядок аргументов. Именно нелинейность и разнообразие вариаций делают изображение органическим, «пламенным», а не полигональным.
Ниже — несколько примеров таких вариаций (обозначим $r=\sqrt{x^2+y^2}$, $\theta=\operatorname{atan2}(x,y)$):
Linear ($V_0$): $V(x,y) = (x, y)$
Sinusoidal ($V_1$): $V(x,y) = (\sin x,\ \sin y)$
Spherical ($V_2$): $V(x,y) = \dfrac{1}{r^2}(x, y)$
Swirl ($V_3$): $V(x,y) = (x\sin(r^2) - y\cos(r^2),\ x\cos(r^2)+y\sin(r^2))$
Horseshoe ($V_4$): $V(x,y) = \dfrac{1}{r}\big((x-y)(x+y),\ 2xy\big)$
Polar ($V_5$): $V = \left(\dfrac{\theta}{\pi},\ r-1\right)$
Handkerchief ($V_6$): $V = \big(r\sin(\theta+r),\ r\cos(\theta-r)\big)$
Heart ($V_7$): $V = \big(r\sin(\theta r),\ -r\cos(\theta r)\big)$
Spiral ($V_9$): $V = \dfrac{1}{r}\big(\cos\theta+\sin r,\ \sin\theta-\cos r\big)$
Julia ($V_{13}$): $V = \sqrt{r}\left(\cos\left(\dfrac{\theta}{2}+\omega\right),\ \sin\left(\dfrac{\theta}{2}+\omega\right)\right),\ \omega \in \{0,\pi\}$ случайно
Всего в классическом наборе Apophysis около 50 базовых вариаций; в разобранной реализации используется 49, включая параметрические — вариации со случайными коэффициентами, которые генерируются один раз на каждый новый узор (перед запуском chaos game), а не на каждую итерацию. Пример — curl:
$$V(x,y) = \frac{1}{t_1^2+t_2^2}\big(x t_1 + y t_2,\ y t_1 - x t_2\big), \qquad t_1 = 1 + p_1 x + p_2(x^2-y^2), \qquad t_2 = p_1 y + 2 p_2 xy$$
где $p_1, p_2$ — случайные параметры, фиксированные для конкретного изображения. За счёт таких параметров одна и та же формула вариации порождает разные семейства форм от запуска к запуску.
В общем виде функция преобразования может использовать сразу несколько вариаций с весовым вектором $v_{ij}$ (ненулевым для нескольких вариаций из всех 49):
$$F_i(x,y) = \sum_j v_{ij}\, V_j(a_ix+b_iy+c_i, d_ix+e_iy+f_i)$$
На практике почти всегда используют один-два ненулевых веса на преобразование — так проще подбирать эстетичные параметры вручную. Рассматриваемая реализация как раз жёстко ограничивается этим частным случаем — не более двух вариаций на преобразование, с весом смешивания $m_i \in (0,1)$:
$$F_i(x,y) = (1-m_i)\, V^{(1)}(x,y) + m_i\, V^{(2)}(x,y)$$
Это увеличивает разнообразие итогового узора без увеличения числа самих аффинных преобразований.
Стоит знать (хотя это не реализовано в разбираемом коде), что в статье Дрейвса функция преобразования расширяется ещё двумя способами. Во-первых, добавляется пост-трансформ $P_i$ — ещё одно аффинное преобразование, применяемое после вариаций: $F_i(x,y) = P_i\big(\sum_j v_{ij}V_j(\dots)\big)$ — он даёт возможность независимо менять «систему координат» каждой вариации. Во-вторых, вводится единый финальный трансформ, применяемый к точке после выбора и применения $F_i$, но уже вне игры в хаос — то есть он не участвует в самой динамике (следующая итерация всё равно берёт координаты до его применения), а действует как «объектив камеры», преобразующий то, что будет напечатано на холст. С учётом пост-трансформов полная параметризация одной функции требует уже 93 параметра, а не 6 — поэтому в редакторах вроде Apophysis у каждого преобразования есть отдельные вкладки для аффинной части, вариаций и пост-трансформа.

сетка 2×3 — одно и то же множество аффинных преобразований, к которым применены разные вариации (linear, spherical, swirl, horseshoe, polar, julia). Наглядно показывает, что именно вариация определяет «стиль» узора при фиксированной геометрии.

На каждой итерации нужно выбрать одно из $n$ преобразований $F_1,\dots,F_n$, причём не равновероятно — у каждого есть вес $w_i > 0$ (в реализации — случайное число из $[0.3,\ 1.0]$), определяющий, насколько часто это преобразование «участвует» в построении узора, а значит — насколько плотным будет соответствующий ему фрагмент изображения. (Некоторые ранние реализации привязывали вес к коэффициенту сжатия каждого преобразования — это ускоряет заполнение при однобитном рендере, — но для нелинейных вариаций коэффициент сжатия не постоянен, и статья прямо советует так не делать; в современных реализациях, включая разбираемую, вес — независимый свободный параметр.)
Стандартный способ выбрать индекс с вероятностью, пропорциональной весу — метод рулетки (roulette wheel selection):
$$p_i = \frac{w_i}{\sum_{j=1}^n w_j}$$

1. Посчитать $W = \sum_j w_j$.
2. Сгенерировать $r \sim \text{Uniform}(0, W)$.
3. Идти по преобразованиям, вычитая $w_i$ из $r$; как только $r$ становится меньше очередного $w_i$ (эквивалентно — накопленная сумма весов превысила $r$), выбрать это преобразование.

У Дрейвса симметрия строится иначе, чем в разбираемом коде, и об этом стоит сказать отдельно. В оригинальном алгоритме симметрию встраивают прямо в систему: в набор преобразований добавляют одну или несколько дополнительных функций, каждая из которых — чистый поворот вокруг начала координат, и эти функции участвуют в игре в хаос наравне со всеми остальными (то есть могут быть выбраны на очередном шаге, и их результат сам становится входом для следующей итерации). Для $k$-кратной поворотной симметрии нужно $k-1$ таких дополнительных функций; при этом веса подобраны не произвольно — каждой из добавленных функций-поворотов даётся вес, равный сумме весов всех остальных функций системы, иначе плотность разных лучей симметрии окажется разной. (Автор поясняет это на примере 3-кратной симметрии: чтобы попасть в «третий» луч, нужно два раза подряд выбрать поворот на 120° — что вдвое менее вероятно, чем один раз выбрать поворот сразу на 240°; поэтому нужны обе функции-поворота, а не одна.) Двусторонняя (диэдральная) симметрия добавляется похоже — функцией, отражающей $x$-координату. Отдельный тонкий момент: при выборе функции-поворота или отражения нельзя обновлять цветовую координату $u$ по общему правилу — иначе после нескольких применений точка возвращается почти в то же место с другим цветом, и палитра «вымывается»; поэтому для симметрийных функций цвет оставляют без изменений.
Реализация, разобранная в этой лекции, устроена проще и дешевле: она не трогает саму игру в хаос и не добавляет никаких новых функций в систему — вместо этого уже посчитанная на шаге итерации точка $(x,y)$ просто откладывается на холст $k$ раз, в повёрнутых на углы $2\pi s/k$ позициях:
$$R_s(x,y) = \big(x\cos\alpha_s - y\sin\alpha_s,\ x\sin\alpha_s + y\cos\alpha_s\big), \qquad \alpha_s = \frac{2\pi s}{k},\ s = 0,\dots,k-1$$
Иначе говоря, здесь симметрия — свойство финального растра, а не самого аттрактора: динамика точки $(x,y)$ по-прежнему асимметрична, симметричны только $k$ её «эхо»-проекций на холсте. У этого подхода два практических следствия по сравнению с методом Дрейвса: во-первых, не нужно вручную балансировать веса — проблема неравной плотности лучей здесь не возникает, потому что каждая проекция получает ровно одну точку за итерацию; во-вторых, цвет точки автоматически одинаков у всех $k$ копий (это буквально одна и та же точка с одним и тем же $u$), так что и проблема «вымывания» цвета сама по себе не возникает. Ценой является то, что получившийся узор — это симметризация одной конкретной (асимметричной) траектории, а не аттрактор действительно расширенной, симметричной системы функций; на глаз разница почти не видна, но с точки зрения теории это два разных объекта.
Каждому преобразованию $F_i$ сопоставлен фиксированный оттенок (hue) $h_i \in [0,1)$, выбираемый случайно один раз при генерации узора. Цвет текущей точки представляет собой экспоненциальное скользящее среднее оттенков всех преобразований, через которые прошла точка:
$$u_{k+1} = \frac{u_k + h_{i_k}}{2}$$
где $i_k$ — индекс преобразования, выбранного на шаге $k$. Так соседние по траектории точки получают близкие, плавно меняющиеся цвета — отсюда характерные цветовые «разводы» на изображениях fractal flame, а не резкие цветовые границы. В общем случае Дрейвс определяет цвет как произвольную функцию-палитру из $[0,1]$ в $(r,g,b)$ (классически — таблица из 256 троек байт); в разбираемой реализации эта палитра зафиксирована как HSV с $S=0.7,\ V=1.0$ и переменным $H=u_k$, затем переводится в RGB по стандартной формуле HSV→RGB.
Ключевое отличие от классической игры в хаос: вместо того, чтобы просто закрашивать пиксель под точкой, для каждого пикселя ведутся четыре независимых аккумулятора:
$$\text{hist}_R[p] \mathrel{+}= r_k, \quad \text{hist}_G[p] \mathrel{+}= g_k, \quad \text{hist}_B[p] \mathrel{+}= b_k, \quad \text{hist}_N[p] \mathrel{+}= 1$$
где $p$ — индекс пикселя, в который спроецировалась точка $x_k$ (с учётом поворотов симметрии из предыдущего раздела), а $(r_k,g_k,b_k)$ — её цвет. В терминологии статьи Дрейвса именно $\text{hist}_N$ и называется альфа-каналом (α) — четвёртым каналом, который вместе с логарифмом даёт итоговую формулу тонирования ниже. $\text{hist}_N[p]$ — это, по сути, оценка плотности вероятности посещения точки $p$ хаотической траекторией.
Точка $x_k$ приближается к аттрактору постепенно; первые несколько итераций она всё ещё находится под влиянием случайно выбранной начальной точки $x_0$, которая почти наверняка не лежит на аттракторе. Решение — отбросить первые $k_0$ итераций и начинать накопление гистограммы только при $k \ge k_0$ (в реализации $k_0 = 50$). Это стандартный приём и в MCMC-семплировании, где его называют burn-in по той же причине.
Не все случайно сгенерированные комбинации аффинных коэффициентов и вариаций дают ограниченную траекторию — вариации вроде tangent ($V=(\sin x/\cos y,\ \tan y)$) или secant содержат особые точки, где значение стремится к бесконечности. Если очередная точка «улетает», практическое решение — перезапустить точку из новой случайной позиции и продолжить игру, не считая эту итерацию частью гистограммы.
Естественная идея — нормировать яркость пикселя линейно по числу попаданий: $\alpha[p] = \text{hist}_N[p] / N_{\max}$. Проблема в том, что распределение плотности в fractal flame имеет чрезвычайно широкий динамический диапазон: в отдельных точках (обычно вблизи неподвижных точек преобразований) плотность на много порядков выше, чем в среднем по изображению. При линейной нормировке почти вся картинка становится либо чёрной, либо перенасыщенной.
Решение Дрейвса — нормировать не саму плотность, а её логарифм:
$$\alpha[p] = \frac{\ln\big(\text{hist}_N[p]+1\big)}{\ln\big(N_{\max}+1\big)}$$
Логарифм «сжимает» большие значения гораздо сильнее, чем маленькие, поэтому в итоговом изображении остаются видны одновременно и плотные, и разреженные области — эффект, схожий с HDR-тональной компрессией в фотографии.

одно и то же накопленное изображение, отрендеренное тремя способами — линейно, логарифмически без гаммы, логарифмически с подобранной гаммой.

Дополнительно применяется степенная гамма-коррекция для точной настройки контраста:
$$\alpha_\gamma[p] = \alpha[p]^{1/\gamma}$$
При $\gamma > 1$ полутона становятся ярче — проявляются более тонкие детали периферии узора; при $\gamma < 1$, наоборот, контраст усиливается, а слабые детали «тонут» в фоне. Автор отмечает, что для некоторых узоров имеет смысл пробовать гамму вплоть до 4 — она дополнительно проявляет мелкие детали, но одновременно усиливает видимый шум, поэтому для больших значений гаммы обычно нужно больше итераций, чтобы шум не был заметен.
Есть и более тонкая проблема: если применять гамма-коррекцию к каждому цветовому каналу независимо, при большой гамме насыщенные цвета «вымываются» в блёклые, пастельные — потому что коррекция подтягивает маленькие значения к единице, стирая разницу между каналами. У Дрейвса за это отвечает отдельный параметр, vibrancy: вместо поканальной коррекции можно посчитать единый масштабный коэффициент по альфа-каналу и применить его сразу ко всем трём каналам, сохранив исходное соотношение между ними (и, соответственно, насыщенность цвета). В разбираемой реализации гамма применяется именно так — общий множитель $\alpha_\gamma[p]$ одинаков для R, G и B:
$$\text{color}_{R}[p] = \min\left(255,\ \frac{\text{hist}_R[p]}{\text{hist}_N[p]} \cdot \alpha_\gamma[p]\right)$$
и аналогично для $G$, $B$ — то есть эта реализация фактически всегда работает в режиме «vibrancy = 1», и проблема выцветания красок в ней не возникает.

Ссылка на оригинальную работу Скотта Дрейвса "The Fractal Flame Algorithm" опубликованную в 2003 году: https://flam3.com/flame_draves.pdf