Фракталы — потрясающая штука. Они всплывают везде в природе, по сути их можно считать в буквальном смысле природной красотой. Самым популярным фракталом является множество Мандельброта — вернее его графическое отображение.
Недавно мне любопытства ради захотелось разобраться, как можно рисовать множество Мандельброта в различных вариациях, а заодно попрактиковаться в работе с языком программирования D. Конечно, рисовать фракталы — не самая сложная задача сама по себе, но от этого не убавился интерес от исследования темы и создания мелких инструментов по пути.
В процессе изучения вопроса порой приходилось закапываться в несколько разных источников, так что я решил совместить приятное с полезным. В этом посте я расскажу всё, что узнал в ходе экспериментов.
 |
| Кадр из Теории Большого Взрыва |
Кстати, этот кадр из The Big Bang Theory использует видео "Eye of The Universe", сгенерированное и загруженное на канал Maths Town. 17 миллионов итераций и зум на 3.4e1091. По этой ссылке есть как само видео, так и сопутствующие материалы и ряд крутых кадров оттуда: https://www.maths.town/videos/eye-of-the-universe-video
Код с инструкциями к выполнению и сборке: https://github.com/leamare/mandelbrot-experiment-dlang
Ну а если вас интересуют красивые картинки, то определённо стоит обратить внимание на галерею в конце!
Началось всё с видео на канале Mathologer, в котором было рассказано про "Буддаброт".
Также сюда можно приплести и другие интересные материалы, в частности видео Veritasium про логистическое отображение или видео про фракталы от 3blue1brown. Но я отвлёкся.
Ну начнём с того, что такое вообще фракталы. Фрактал — это, грубо говоря, множество, которое отдельных своих регионах похоже само на себя. Очень частый пример — деревья. Если взять кусочек дерева — ветку — и посмотреть отдельно, то она будет похожа на маленькое дерево со своими маленькими ветками, и продолжать так можно очень долго.
Фракталов существует много (как и определений самому понятию). Но самым известным можно считать множество Мандельброта. Рассчитывается всё довольно просто по формуле
где начальным значением z будет 0, а c — комплексное число. Если при повторении этой операции много раз результат будет стремиться к бесконечности, значит число не является частью множества. Если же какой-то лимит значения существует — число является частью множества. Вот так вот просто.
Ну и ещё есть верхний лимит количества итераций ("dwell") — количество раз, которое будут производиться операции над точкой. Чем больше итераций, тем точнее результат и тем больше деталей на графике.
Есть ещё фрактал Жюлиа, и разные его состояния можно найти на границе множества Мандельброта, как в каталоге, но это не так важно сейчас.
Но если очень интересно, есть очень прикольное видео, наглядно показывающее связь множеств Мандельброта и Жюлеа: https://www.youtube.com/watch?v=dekA3GNm6Qk
Есть разные вариации этого множества. Самые известные — "корабль" (в нём в степень возводится модуль комплексного числа) и "мультиброт" (где стоит любая степень помимо 2). Все техники, впрочем, применимы и к ним ровно в той же степени, так что пока что я буду говорить только об основном.
 |
| Семейство Бротов: Мандельброт, Мультиброт (степень 4), Мультиброт (степень -2) и Корабль |
Определить какую-то зону множества можно по сути прямоугольником (координаты углов), либо точкой центра с "радиусом" камеры.
Существует несколько различных методик, которыми можно нарисовать множество Мандельброта (ну и фрактал Жюлиа тоже, да и вообще всё). Можно просто нарисовать чёрным точки, попадающие в множество, а белым — не попадающие. Можно их покрасить разными цветами в зависимости от того, как много итераций нужно сделать, чтобы достигнуть "лимита". Можно ещё определять цвета точек на основе того, как часто путь на "выход" пролегает через эту точку. Можно даже построить трёхмерную версию (или четырёхмерную, но такие вещи лучше поищите на ютубе).
Я обращал внимание на три методики: основная (цвет на основе количества итераций), "буддаброт" и "анти буддаброт" (о них позже).
Покрас по итерациям
Ну начать стоит с того, что количество итераций целочисленно, и графики на основе целочисленного значения получаются... ну не самые красивые, скажем так. Хочется, чтобы цвета сменялись плавно.
 |
| Результат самого простого подхода к покрасу в самой первой версии |
iter_d = iterations + 1 - (ln(2) / |Z[n]|) / ln(2)Ну или вариант в виде кода:
const double log_zn = log(Zr*Zr + Zi*Zi) * 0.5;
const double nu = log( log_zn / log(2.0) ) / log(2.0);
iter_d = 1 + to!double(iter) - nu;Но, само собой, можно немного ускорить код, заранее посчитав log(2.0), но это уже другая история.
Но что делать с цветами?
Можно подойти к вопросу с двух сторон: использовать заранее заданную палитру-градиент или считать цвет целиком.
В первом случае подход простой: считаем, между какими цветами точка должна находиться, после чего берём соответствующее промежуточное значение.
 |
| Градиент ultrafractal, X=-1.9990956182679107, Y=2.9177492662331486e-09, R=1.304734098539484e-12, dwell=5000, размер палитры 250 |
Есть некоторые популярные градиенты, которыми пользуются крутые ребята, вроде градиента "Earth and Sky" (он же ultrafractal, используется на картинках с википедии), если нужны ключевые точки этих градиентов — можно заглянуть в мой код. Ключевые точки у меня захардкожены прямо в mandel.d, там же есть несколько особых градиентов, которые я добавил сам.
Второй случай интереснее. Тут тоже можно решить проблему разными способами. Например, можно просто "назначить" каждой точке соответствующий "серый" цвет (RGB цвет с одинаковыми значениями на всех каналах). При этом, чтобы выглядело красивенько, нужно бы ещё чтобы градиент шёл не только по нарастающей, но и по убывающей.
// например так
auto v = 2*(iter_d/paletteSize) % 2; // исходное значение, которое пригодится при рассчётах
const auto vc = v > 1 ? 2-v : v; // "повёрнутый" цвет |
| Grayscale |
Если же хочется, чтобы всё было окрашено одним цветом, можно просто ставить в канал "основного" цвета полное значение, а в остальных каналах играться с множителями.
У меня в коде (в том же модуле mandel, функция pixelcolor()) можно посмотреть, какие я ставил множители для Ocean Blue / Fire Red, если очень хочется.
 |
| Ocean Blue |
Самым интересным становится, пожалуй, окрас через HSV цвета. Если в RGB каждый цвет делится на красный, зелёный и синий каналы, то в HSV первое значение определяет тон (и обычно вообще выражается как "градус" или положение на окружности), второе — насыщенность, а третье — непосредственно значение (условно говоря "яркость"). Но, само собой, просто скормить значение v прямо сюда было бы глупо, понадобятся дополнительные модификации, зависящие ещё и от реализации библиотеки, через которую происходит работа с HSV цветами. У меня получилось достичь более-менее приятных глазу результатов с такими манипуляциями:
const auto vc = v > 1 ? 2-v : v;
auto vl = 0.25+vc*2;
auto vl = 0.25+vc*2;
auto vs = 0.75+vc*2;
auto c = hsv(
360.0*v,
vs > 1 ? 1 : vs,
vl > 1 ? 1 : vl
);  |
| HSV |
Ну и ещё одна проблема, которую стоит учесть — сглаживание. Чаще всего для сглаживания просто считают изображение в разрешении в два раза большем, чем нужно, после чего понижают разрешение вдвое, усредняя цвета соседних точек.
Если вы искали какие-нибудь клёвые градиенты (или "тот самый градиент с википедии"), то у меня в коде есть несколько таких (прямо "вшиты" в mandel.d):
- ultrafractal — он же "Earth and Sky"
- seashore
- fire и oceanid — аналоги red и blue в виде заранее определённых градиентов, а не полностью вычисляемые
- cnfsso — градиент, который мне подсказала Настюша
- acid — ещё один градиент от Настюши, кислота по-дизайнерски
- softhours — и ещё один (хотя тут вышли слишком пастельные цвета, плохо видно)
Ну и есть также реализации для HSV покраса, ЧБ варианта и Синей/Красной палитры.
Буддаброт
Вообще говоря есть два варианта этого графика: включающий точки внутри самого множества и исключающий их. Первый называют "Анти Буддабротом", второй — "Буддабротом".
Само название происходит из того, что график множества, будучи повёрнутым на 90 градусов, похож на сидящего Будду.
 |
| Будда и Антибудда |
Впервые Буддаброт нашла Мелинда Грин, и это поистине завораживающая картина.
По этой ссылке можно почитать оригинальный текст по теме, а по этой ссылке будет статья в википедии
Как я уже писал выше, считаются значения точек Буддаброта путём отслеживания "маршрутов" точек во время итераций. Методов рассчёта есть несколько. Одним из самых распространённых является "случайный бросок", когда прослеживаются маршруты случайных точек на плоскости.
Я же использовал более очевидный вариант: при отрисовке основного графика множества всё равно считаются значения для всех точек, так почему бы не использовать именно эти просчёты, чтобы "заполнить" данные о буддаброте?
Есть ещё один вариант: взять изображение с множеством Мандельброта, после чего выделить края множества, а далее итерировать только точки внутри этих краёв. Такая оптимизация сильно ускоряет вычисления, при условии, что у вас не стоит целью работать с каждой точкой и без этого. В моём случае такое ускорение было бы несколько бессмысленно, так как я как раз считаю одновременно оба графика.
Однако вычисление Буддаброта всё же занимает больше времени, а главная причина тому — постоянные обращения к общему хранилищу данных о "проитерированных" точках, что тормозит многопоточность.
Это мой косяк, конечно, но мне было лень сильно задумываться на этот счёт и искать иной подход.
Но есть тут один нюанс: результирующее изображение неизбежно будет очень шумным. Бороться с этим можно несколькими путями:
- добавить больше точек для вычислений (сделать дополнительный "случайный бросок" или добавить точек вне основного вычисляемого множества)
- сделать очень большое изображение и понижать разрешение (это мой вариант)
- использовать алгоритмы подавления шума (алгоритмов много разных, и я для них слишком глупый, но некоторые люди подходили к вопросу проще и брали, например, OptiX от NVIDIA и применяли для подавления шума)
Я использую комбинацию из последних двух пунктов. Изображения с буддабротом, которые есть в этом тексте сделаны именно так.
 |
| Тестовый Небулаброт с маленьким разрешением и исходным уровнем шума |
Ах, да, туманности. Самая интересная часть — как можно покрасить Буддаброт? Ведь по определению он будет чёрно-белым. И на самом деле уже в этом месте есть небольшая проблема: во-первых, яркость разных итераций будет отличаться, а с ростом максимального количества итераций график вообще становится тяжелоразличимым. Несколько нормализовать яркость можно разными способами, я решил для этого брать корень четвёртой степени от отношения между итоговым количеством итераций и максимальным, что не решает проблему полностью, но делает картинку достаточно яркой, чтобы детали можно было различить.
 |
| Ещё один тестовый небулаброт, показывающий уровни яркости без нормализации (не обращайте внимание на визуальные артефакты) |
При этом даже после таких манипуляций всё равно придётся отдельно поиграться, чтобы сделать туманность. А что за туманность вообще?
По сути это просто очередное отображение того же самого графика. За основу берутся три Буддаброта с разным максимальным количеством итераций (желательно кратные, например 50, 500 и 5000 или 100, 1000 и 10000). Далее изображения идут в синий канал (с самым маленьким количеством максимальных итераций), в зелёный (средним) и красный (максимальным). Итоговое изображение в целом и так выглядит красиво, но чтобы всё было совсем уже классно, придётся подкрутить параметры цветовых каналов вручную.
Забавный факт: этот метод взят из космической фотографии, где вместо количества итераций используются разные уровни выдержки.
 |
| Буддаброты, полученные с использованием разного количества итераций |
 |
| Основное множество, 100/1000/10000 итераций |
Чтобы получить самые приятные глазу буддаброты, я использовал базовое число итераций 50 и разрешение 32k x 32k пикселей, уменьшенное до 8k x 8k (хотя для некоторых я использовал 16k x 16k).
 |
Шея Будды (X=-1.15 Y=0 R=0.125), немного окрашена в синий, 50/500/5000 итераций |
 |
| Малый остров (X=-1.25275 Y=-0.343 R=0.0025), схожий с тем, что был на странице Мелинды, 500/5000/50000 итераций |
Язык D и другие технические штучки
Мне вообще всегда нравился C, но не очень нравилась возня с указателями и всем вот этим вот. Когда я в первый раз попробовал D несколько лет назад (точнее говоря, лет 10 уже), я был в восторге от того, насколько приятнее работать на D. Я попробовал экспериментировать с языком ещё тогда, но делать с ним было особо нечего.
С тех пор язык успел преобразиться, обрасти сообществом и функционалом, и в целом опыт написания программ на D напоминает смесь C, Python, Go и JavaScript — каждый по своим причинам, но главное — всё относительно просто и понятно, при этом язык хорошо себя чувствует на "низком уровне".
Но с непопулярностью приходят и некоторые неудобства. Например, отсутствие каких-либо библиотек (которые просто ещё не сделали), либо сломанность других (которые просто не поддерживали пять лет). У некоторых библиотек нет адекватной документации, а за описанием некоторых вещей приходится и вовсе лезть в документацию компилятора. Но это обычно происходит редко, в особо специфичных случаях, и только добавляет интерес.
На самом деле эти проблемы есть у всех языков, но у более популярных до такого обычно редко доходит.
Что же касается реализации, тут было несколько пунктов, про которые я думал и с которыми столкнулся, и про каждый из них я бы хотел поговорить.
Параллелизм
Это потрясающе. Буквально ОДНО изменение в коде в ОДНОЙ строчке и все вычисления равномерно делятся между ядрами. Восхитительно.
Ну, конечно, тут есть свои нюансы. Например, если постоянно обращаться к каким-то внешним данным, толку от параллельных вычислений не будет совсем.
Live renderer (SFML)
На самом деле это не первый мой такой "эксперимент" на D, и пара прошлых уже использовали в своей основе SFML для рендеринга графики на экране.
Но, во-первых, эта библиотека ориентирована на C++, а в D работает через биндинги для C, которые сделаны энтузиастами (и слабо документированы).
Во-вторых, в какой-то момент я решил, что рендеринг картинки занимает достаточно долгое время (вернее не рендеринг, а просчёт точек), и а для ускорения с GPU я не имел опыт общения с шейдерами (и имею некоторые сомнения, что это сильно ускорит процесс по ряду причин). К тому же, мне хотелось рисовать возмутительно большую картинку, а это сильнее замедляет процесс вычислений.
Dlib
Так что я пошёл по иному пути. Я начал рендерить картинку в файл, а верным помощником в этом деле стала библиотека Dlib от Gecko0307 — человека из сообщества D-разработчиков, а также разработчика игр (и делающего свой игровой движок на D).
https://github.com/gecko0307/dlib/
Формат BMP
Кстати про изображения. Изначально я хотел сохранять всё "в чистом BMP", без сжатия в каком-либо виде. Проблема, с которой я столкнулся — если разрешение картинки слишком большое, "количество" точек просто не влезает в заголовок файла, потому что формат не был готов к таким извращениям, так что в итоге он становится нечитабельным.
Комплексные числа
Вообще говоря в языке D (как и везде уже наверное) есть встроенная библиотека для работы с комплексными числами (и даже встроенные типы для этого), но идти по такому пути было бы слишком легко, так что все вычисления с возведением комплексных чисел в степень я писал сам. Ну а больше мне ничего и не нужно было.
Но вообще в коде у меня есть кусок, который использует готовую реализацию (он закомментирован, использовался для сверки результатов).
Числа высокой точности
Вообще говоря тип real, который имеется в языке, зависит от максимального размера памяти под число с плавающей точкой, поддерживаемого процессором на системе. И, как водится, из-за того, как типы данных с плавающей точкой работают, есть некоторые ограничения, после которых точность начинает страдать, и ниже уже никак не опуститься.
В разных реализациях я натыкался на похожий лимит (и всегда он был примерно в районе e-14 и при большом количестве перечисленных знаков после запятой). Если перейти эту границу, генерируемые изображения просто будут либо генерироваться безумно долго, либо покрываться артефактами из-за проблем "округления".
Решить проблему можно было бы используя реализацию "decimal" типа, и можно было бы даже самому её написать, НО
- вычисления с decimal занимают ЗНАЧИТЕЛЬНО больше времени
- мне лень
- есть несколько готовых реализаций, но все они заброшены, а единственная более-менее функциональная только-только начинает переезд в группу репозиториев dlang-community (и она тоже не работает на данный момент)
Так что я забил на высокую точность, хотя мне конечно хотелось сделать крутую часовую анимацию с погружением на e-1024, как в этом видео от Maths Town, которым я вдохновлялся.
Долгие вычисления и сериализация прогресса
В какой-то момент я понял, что вычисление огромных буддабротов занимает очень много времени, так что пришла идея как-то сохранять прогресс. Но "повесить обработчик" на SIGTERM/SIGKILL или любой другой вариант прекращения работы не получится, надо искать иной вариант.
Ну для начала надо найти способ относительно быстро и "дёшево" сериализовать данные в бинарном виде. Изначально я хотел использовать какую-нибудь из реализаций msgpack, но они (естественно) не очень хорошо работали с тем, что я хотел сделать. Так что выбор пал на библиотеку cerealed.
Ну а в остальном реализация получилась несколько топорная, но подходящая для моих нужд: выделяется размер блока, блок выполняется параллельно, как и положено, после чего результаты сохраняются на диск, а после сохранения начинается следующий блок. Это позволяет избежать коллизий и получить хоть какой-то толчок при старте пересчёта.
Такой подход, конечно, сильно затормаживает выполнение, потому что сохранение занимает время, но с другой стороны это нужно в первую очередь для ОЧЕНЬ долгих задач, которые иногда хочется прервать (чтобы заняться чем-то ещё), а потом продолжить, так что такие потери допустимы.
Анимации, рендеринг по кускам и очередь
Тут особо нечего сказать. Мне в какой-то момент захотелось генерировать прикольные анимации, а заодно ещё генерировать изображения "чанками".
Базовые механизмы деления изображения на части (или вернее сказать подготовки зон и центральных точек этих чанков) или анимаций (не считая того факта, что итерации, размер палитры и радиус должны расти используя экспоненту и разницу логарифмов, линейное увеличение тут не поможет) сами по себе простые.
Интереснее было другое: чтобы всё это работало элегантно и без лишних костылей, я переделал изначальную версию рендера, сделав для него очередь задач, в которую по сути и закидываются все запросы на изображения: будь то обычный набор параметров, кадры анимации или секции большого изображения.
Ах, да, ещё я для этого сделал чтение описания задач из JSON файла, в котором указываются все-все-все нужные параметры (координаты точек, радиус камеры, тип фрактала и что только захочешь).
Разные наблюдения и способы улучшить код
Ну о том, что небулаброты в любом случае приходилось подстраивать вручную я сказал. Конечно, можно было бы сделать автоматическое комбинирование каналов и поиграться дальше с нормализацией цвета и всё такое, но мне было лень.
Нормализация яркости, кстати, делает очень сложным использование "генерации по частям" для буддабротов из-за разных уровней яркости у каждого из кусков (придётся подстраивать руками при склейке).
Изначально я пару раз налажал с округлением при конвертации координат, что привело к очень смешным артефактам (потому что я спросонья написал floor вместо round).
Для поиска интересных координат, которые можно отрисовать, я бродил по сети. Много интересных точек я нашёл тут: http://math.hws.edu/eck/js/mandelbrot/MB-info.html
Другие же писал сам.
Точность страдает при глубоком погружении, и мне бы хотелось это исправить, чтобы погрузиться до e-1024. Возможно в будущем вернусь к этой идее и допилю.
Исправить проблему можно было бы использованием чисел высокой точности, как я писал раньше, можно даже реализовать тип Decimal из стандарта IEEE-754-2008.
В целом изобрести какой угодно костыль вокруг этой проблемы не будет сложным, НО это займёт время, а мне его не очень хочется тратить на это. Из готовых реализаций мне на глаза попались несколько decimal библиотек для D, но они были либо не до конца рабочими (например не поддерживали деление и умножение), либо сломанными (из-за отсутствия поддержки). При этом у ряда из них также были описания с бенчмарками, которые показывали значительное проседание по скорости вычислений (до 20 раз).
  |
| Пара примеров того, как могут выглядеть артефакты, полученные из-за округления |
Рассчитать автоматически необходимое число итераций для конкретной точки практически невозможно. Есть либо примерные приближения, которые быстро становятся бесполезными, либо есть методы, для которых надо посчитать с максимальным количеством итераций ближайшие точки в окрестности, чтобы потом решить, а какое количество итераций считать за приемлемое. Есть и другие методы, основанные на этой идее, но мне они не понравились, так что я решил забить на эту идею.
При генерации видео, как я уже писал выше, радиус камеры и количество итераций (ну и палитру) надо считать через разницу логарифмов начального и конечного значений (ну а потом прибавлять полученную дельту к логарифму начальной точки и считать экспоненту). Линейно тут не посчитаешь.
Рассчёт негативных степеней комплексных чисел оказался забавной задачей: положительные степени можно считать без проблем тригонометрически, а для отрицательных нужно сначала посчитать положительную степень, а потом "поделить один на результат" — по сути поделить каждую из составляющих на модуль (Zr² + Zi²) результирующего комплексного числа.
Ещё одна забавная деталь: чтобы получить адекватную картинку мультиброта с негативной степенью, стоит "обратить" число итераций — вычесть из максимального числа итераций результирующее число.
Мне хотелось сделать ещё автоматический поиск значений из множества Жюлеа, но стало лениво. В этой статье есть немного информации по теме (и ещё крутые картиночки):
https://nextjournal.com/lazarus/aggregating-values-to-the-mandelbrot-and-julia-sets
Есть ещё интересные статьи про визуализации множества Мандельброта (и Буддаброты), а также с описанием разных оптимизаций, позволяющих ускорить процесс. Самой интересной статьёй можно считать этот материал от Дэвида Араманта о том, как он считал 68.7 гигапиксельный буддаброт: https://davidaramant.github.io/post/the-buddhabrot-part-5/
Там же есть общая информация, как это всё считать и рисовать, а также про нормализацию яркости через tone mapping и разные оптимизации.
Также про покрас Буддаброта и оптимизации: https://mathematica.stackexchange.com/questions/89458/how-to-make-a-nebulabrot
Шустрый генератор небулабротов на основе WebGL и шейдеров: https://github.com/donghaoren/buddhabrot-renderer (можно заметить, что изображение становится "лучше" постепенно, а с погружением в множество генерация изображения сильно замедляется; но зато тут можно поиграться с разными переменными и даже вращать множество)
А вот в этой статье предлагается иной подход — использование мутаций и алгоритма Метрополиса Гастингса: http://www.steckles.com/buddha/
Галерея
Ну и куда же без галереи картиночек! И даже два видео!
В конце я оставляю несколько самых красивых результатов, два видео и ссылку на папку в Google Drive, где можно посмотреть на все эти картиночки в исходном качестве. Ну а если есть желание воспроизвести результаты, в именах файлов можно заметить значения всех необходимых переменных.
Может быть однажды и мои картиночки на википедии будут висеть, хехе
