Шёл уже последний час этого воскресенья, я уже думал идти спать, но добрый sourcerer прислал мне картинку с моего заброшенного сайта, которую можно увидеть ниже, и текст «красиво!». Эти картинки я

Ночь фракталов

Математика, Алгоритмы*

Шёл уже последний час этого воскресенья, я уже думал идти спать, но добрый sourcerer прислал мне картинку с моего заброшенного сайта, которую можно увидеть ниже, и текст «красиво!». Эти картинки я рисовал лет пять назад, с помощью т. н. алгоритма времени убегания, но для применимости данного алгоритма, нужно уметь для заданного набора преобразований разбивать плоскость на регионы, тогда я не придумал, как это сделать, и больше к этому алгоритму не возвращался. Но сейчас я сразу сообразил, что делать, и написал Диме: «Сначала Random IFS, потом kNN, а затем Escape-Time Algorithm!»

Под рукой у меня был только старый нетбук, который мне дали друзья на время, пока мой ноутбук в ремонте. Дима мне ещё что-то говорил, я ему что-то отвечал, но у меня уже в голове писался код, и я искал на нетбуке хоть какой-нибудь компилятор или интерпретатор и нашёл C++ Builder 6! После этого я понял, что утро я встречу наедине с борландовским компилятором. Через пять часов я отправил Диме новых картинок, но он, как нормальный человек, давно спал…

Итак, немножко формул. Представим, что у нас есть конечный набор преобразований плоскости Ti: R ² → R ², i = 1,..., k. Для произвольного множества E определим T (E) = T ₁(E) ∪… ∪ T _k(E), т. е. подействуем каждым преобразованием на множество E, а результаты объединим. Можно доказать, что если отображения Ti были сжимающими, то последовательность E, T (E), T (T (E)),… сойдётся к некоторому множеству F для любого непустого компактного E. Данная конструкция и известна как система итерируемых функций.

Например, если в качестве непустого компакта взять смайлик, и рассмотреть три преобразования, каждое из которых является композицией сжатия и сдвига в i -ую вершину правильного треугольника, то первые итерации будут выглядеть так, а в пределе получится треугольник Серпинского:

Обычно вместо прямого вычисления последовательности E, T (E), T (T (E)),… для построения фракталов используют т. н. «игру в хаос», которая заключается в следующем. Выберем произвольную точку z ₀ на плоскости, далее выберем случайно преобразование Ti ₁ и вычислим z ₁= Ti ₁(z ₀), далее снова случайно выберем Ti ₂ и вычислим z ₁= Ti ₂(z ₀), и т. д. Можно показать, что всё будет хорошо, и множество полученных точек будет в некотором смысле приближать множество F, определённое выше. На этот алгоритм я ниже буду ссылаться как на Random IFS.

z = (0, 0)

for (i = 0; i < maxIterNum; ++i) {

cl = random([p1,..., pk]) // pi -- вероятность, с которой выбираем преобразование Ti.

z = T[cl](z)

if (i > skipIterNum) { // Первых несколько итераций могут быть достаточно далеко от аттрактора.

draw(z)

}

}

Теперь самое время перейти к описанию алгоритма времени убегания. Пусть у нас для начала есть одно преобразование плоскости f. Для каждой точки z плоскости начнём вычислять последовательность f (z), f (f (z)), f (f (f (z))),… до тех пор пока либо число итераций не превысит некоторого заданного числа N, либо пока норма числа z не станет больше некоторого числа B. После этого цвет точки выбираем в соответствии с количеством произведённых итераций.

for (y = y1; y < y2; y += dy) {

for (x = x1; x < x2; x += dx) {

z = (x, y);

iter = 0;

while (iter < N && ||z|| < B) {

z = f(z)

iter += 1;

}

drawPixel((x, y), color(iter))

}

}

Если на время представить что наша плоскость является комплексной, а преобразование f (z) равно z ² + c, то в результате работы этого алгоритма мы получим фрактальное множество Жюлиа. Более подробно про это можно прочитать в хабрастатье«Построение множества Жюлиа» хабрапользователя mephistopheies.

Пусть теперь у нас есть система итерируемых функций, заданная набором обратимых сжимающих преобразований плоскости T1,..., Tk. Пусть F — это аттрактор этой системы.

Дополнительно предположим, что множество F можно разбить так, что Ti (F) ∩ Tj (F) = ∅, i!= j (это предположение далеко не всегда выполняется). Разобъём всю плоскость R ² на куски A ₁,...., Ak так, что Ti (F) является подмножеством Ai для всех i. Теперь определим функцию f (z) кусочно: на множестве Ai положим f (z) = Tk ⁻¹(z) для всех i.

Например, для треугольника Серпинского рассмотрим такое разбиение (тут есть небольшие проблемы с тремя точками, но закроем на них глаза).

А теперь самый главный вопрос, что получится, если алгоритм времени убегания применить к построенной таким образом функции f?

Давайте посмотрим:

Получился симпатичный такой треугольник Серпинского!

Оказывается это не случайность. Ещё пару примеров:

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

А именно, сначала с помощью Random IFS генерируем несколько тысяч точек, при этом для каждой точки запоминаем номер преобразования, с помощью которого она была получена. Затем во время работы EscapeTimeAlgorithm для каждого пикселя определяем область, в которую он попадаем с помощью 1NN.

Например, для такой звёздочки 1NN даёт следующее разбиение на четыре куска:

Собирая вместе, получим:

points = RandomIFS(Ts)

classifier = kNN(points);

for (y = y1; y < y2; y += dy) {

for (x = x1; x < x2; x += dx) {

z = (x, y)

iter = 0

while (iter < maxIterNum && ||z|| < sqrdBound) {

cl = classifier.getClass(z);

z = T[cl].applyInvert(z);

iter += 1;

}

draw((x, y), color(iter))

}

}

Ещё несколько картинок.

Вот и всё. Напоследок два замечания.

Во-первых, внимательный читатель возможно задался вопросом, раз фракталы, которые строятся с помощью Random IFS, можно построить с помощью алгоритма времени убегания, то можно ли множество Жюлиа построить с помощью Random IFS? Оказывается можно, нужно просто обратить отображение f (z) = z ² + c, вспомнив, как извлекается корень из комплексного числа. (Правда при применении этого метода для построения изображений множества Жюлиа возникают большие трудности.)

x = z0.re

y = z0.im

for (i = 0; i < N; ++i) {

x -= c.re;

y -= c.im;

len = sqrt(x*x + y*y);

ang = atan2(y, x);

if (rand() % 2 == 0) { // Тут нужно что-нибудь и поинтереснее.

x = sqrt(len) * cos(ang / 2);

y = sqrt(len) * sin(ang / 2);

} else {

x = sqrt(len) * cos(ang / 2 + M_PI);

y = sqrt(len) * sin(ang / 2 + M_PI);

}

draw(x, y)

}

Во-вторых, в статье было рассказано что происходит, если вы хотите узнать почему так происходит, то рекомендую книгу M. Barnsley «Fractals Everywere».

Дата добавления: 2015-10-21