Ночь фракталов
- пятница, 17 октября 2014 г. в 03:11:29
http://habrahabr.ru/post/240417/
Шёл уже последний час этого воскресенья, я уже думал идти спать, но добрый sourcerer прислал мне картинку с моего заброшенного сайта, которую можно увидеть ниже, и текст «красиво!». Эти картинки я рисовал лет пять назад, с помощью т. н. алгоритма времени убегания, но для применимости данного алгоритма, нужно уметь для заданного набора преобразований разбивать плоскость на регионы, тогда я не придумал, как это сделать, и больше к этому алгоритму не возвращался. Но сейчас я сразу сообразил, что делать, и написал Диме: «Сначала Random IFS, потом kNN, а затем Escape-Time Algorithm!»
Под рукой у меня был только старый нетбук, который мне дали друзья на время, пока мой ноутбук в ремонте. Дима мне ещё что-то говорил, я ему что-то отвечал, но у меня уже в голове писался код, и я искал на нетбуке хоть какой-нибудь компилятор или интерпретатор и нашёл C++ Builder 6! После этого я понял, что утро я встречу наедине с борландовским компилятором. Через пять часов я отправил Диме новых картинок, но он, как нормальный человек, давно спал…
Итак, немножко формул. Представим, что у нас есть конечный набор преобразований плоскости Ti: R2 → R2, i = 1, ..., k. Для произвольного множества E определим T(E) = T1(E) ∪… ∪ Tk(E), т. е. подействуем каждым преобразованием на множество E, а результаты объединим. Можно доказать, что если отображения Ti были сжимающими, то последовательность E, T(E), T(T(E)),… сойдётся к некоторому множеству F для любого непустого компактного E. Данная конструкция и известна как система итерируемых функций.
Например, если в качестве непустого компакта взять смайлик, и рассмотреть три преобразования, каждое из которых является композицией сжатия и сдвига в i-ую вершину правильного треугольника, то первые итерации будут выглядеть так, а в пределе получится треугольник Серпинского:
Обычно вместо прямого вычисления последовательности E, T(E), T(T(E)),… для построения фракталов используют т. н. «игру в хаос», которая заключается в следующем. Выберем произвольную точку z0 на плоскости, далее выберем случайно преобразование Ti1 и вычислим z1= Ti1(z0), далее снова случайно выберем Ti2 и вычислим z1= Ti2(z0), и т. д. Можно показать, что всё будет хорошо, и множество полученных точек будет в некотором смысле приближать множество F, определённое выше. На этот алгоритм я ниже буду ссылаться как на Random IFS.
Теперь самое время перейти к описанию алгоритма времени убегания. Пусть у нас для начала есть одно преобразование плоскости f. Для каждой точки z плоскости начнём вычислять последовательность f(z), f(f(z)), f(f(f(z))),… до тех пор пока либо число итераций не превысит некоторого заданного числа N, либо пока норма числа z не станет больше некоторого числа B. После этого цвет точки выбираем в соответствии с количеством произведённых итераций.
Если на время представить что наша плоскость является комплексной, а преобразование f(z) равно z2 + c, то в результате работы этого алгоритма мы получим фрактальное множество Жюлиа. Более подробно про это можно прочитать в хабрастатье «Построение множества Жюлиа» хабрапользователя mephistopheies.
Пусть теперь у нас есть система итерируемых функций, заданная набором обратимых сжимающих преобразований плоскости T1, ..., Tk. Пусть F — это аттрактор этой системы.
Дополнительно предположим, что множество F можно разбить так, что Ti(F) ∩ Tj(F) = ∅, i != j (это предположение далеко не всегда выполняется). Разобъём всю плоскость R2 на куски A1, ...., Ak так, что Ti(F) является подмножеством Ai для всех i. Теперь определим функцию f(z) кусочно: на множестве Ai положим f(z) = Tk−1(z) для всех i.
Например, для треугольника Серпинского рассмотрим такое разбиение (тут есть небольшие проблемы с тремя точками, но закроем на них глаза).
А теперь самый главный вопрос, что получится, если алгоритм времени убегания применить к построенной таким образом функции f?
Давайте посмотрим:
Получился симпатичный такой треугольник Серпинского!
Оказывается это не случайность. Ещё пару примеров:
В этих примерах соответствующее разбиение плоскости не сложно задать аналитически с помощью булевых комбинаций кругов и полуплоскостей, используя метод пристального вглядывания. Но часто простых условий угадать не удаётся. Поэтому вместо угадывания мы научим компьютер определять разбиение самостоятельно. В этом нам поможет метод ближайшего соседа.
А именно, сначала с помощью Random IFS генерируем несколько тысяч точек, при этом для каждой точки запоминаем номер преобразования, с помощью которого она была получена. Затем во время работы EscapeTimeAlgorithm для каждого пикселя определяем область, в которую он попадаем с помощью 1NN.
Например, для такой звёздочки 1NN даёт следующее разбиение на четыре куска:
Собирая вместе, получим:
Ещё несколько картинок.
Вот и всё. Напоследок два замечания.
Во-первых, внимательный читатель возможно задался вопросом, раз фракталы, которые строятся с помощью Random IFS, можно построить с помощью алгоритма времени убегания, то можно ли множество Жюлиа построить с помощью Random IFS? Оказывается можно, нужно просто обратить отображение f(z) = z2 + c, вспомнив, как извлекается корень из комплексного числа. (Правда при применении этого метода для построения изображений множества Жюлиа возникают большие трудности.)
Во-вторых, в статье было рассказано что происходит, если вы хотите узнать почему так происходит, то рекомендую книгу M. Barnsley «Fractals Everywere».
(Мгновения исходного кода можно найти на гитхабе.)
Шёл уже последний час этого воскресенья, я уже думал идти спать, но добрый sourcerer прислал мне картинку с моего заброшенного сайта, которую можно увидеть ниже, и текст «красиво!». Эти картинки я рисовал лет пять назад, с помощью т. н. алгоритма времени убегания, но для применимости данного алгоритма, нужно уметь для заданного набора преобразований разбивать плоскость на регионы, тогда я не придумал, как это сделать, и больше к этому алгоритму не возвращался. Но сейчас я сразу сообразил, что делать, и написал Диме: «Сначала Random IFS, потом kNN, а затем Escape-Time Algorithm!»
Под рукой у меня был только старый нетбук, который мне дали друзья на время, пока мой ноутбук в ремонте. Дима мне ещё что-то говорил, я ему что-то отвечал, но у меня уже в голове писался код, и я искал на нетбуке хоть какой-нибудь компилятор или интерпретатор и нашёл C++ Builder 6! После этого я понял, что утро я встречу наедине с борландовским компилятором. Через пять часов я отправил Диме новых картинок, но он, как нормальный человек, давно спал…
Итак, немножко формул. Представим, что у нас есть конечный набор преобразований плоскости Ti: R2 → R2, i = 1, ..., k. Для произвольного множества E определим T(E) = T1(E) ∪… ∪ Tk(E), т. е. подействуем каждым преобразованием на множество E, а результаты объединим. Можно доказать, что если отображения Ti были сжимающими, то последовательность E, T(E), T(T(E)),… сойдётся к некоторому множеству F для любого непустого компактного E. Данная конструкция и известна как система итерируемых функций.
Например, если в качестве непустого компакта взять смайлик, и рассмотреть три преобразования, каждое из которых является композицией сжатия и сдвига в i-ую вершину правильного треугольника, то первые итерации будут выглядеть так, а в пределе получится треугольник Серпинского:
Обычно вместо прямого вычисления последовательности E, T(E), T(T(E)),… для построения фракталов используют т. н. «игру в хаос», которая заключается в следующем. Выберем произвольную точку z0 на плоскости, далее выберем случайно преобразование Ti1 и вычислим z1= Ti1(z0), далее снова случайно выберем Ti2 и вычислим z1= Ti2(z0), и т. д. Можно показать, что всё будет хорошо, и множество полученных точек будет в некотором смысле приближать множество 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) равно z2 + c, то в результате работы этого алгоритма мы получим фрактальное множество Жюлиа. Более подробно про это можно прочитать в хабрастатье «Построение множества Жюлиа» хабрапользователя mephistopheies.
Пусть теперь у нас есть система итерируемых функций, заданная набором обратимых сжимающих преобразований плоскости T1, ..., Tk. Пусть F — это аттрактор этой системы.
Дополнительно предположим, что множество F можно разбить так, что Ti(F) ∩ Tj(F) = ∅, i != j (это предположение далеко не всегда выполняется). Разобъём всю плоскость R2 на куски A1, ...., Ak так, что Ti(F) является подмножеством Ai для всех i. Теперь определим функцию f(z) кусочно: на множестве Ai положим f(z) = Tk−1(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) = z2 + 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».
(Мгновения исходного кода можно найти на гитхабе.)