Исследуем множество Мандельброта

Изменено: . Разместил: FadeDemon | Меню
Количество просмотров5`855 / FD981-6-0

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


Что же за зверь такой, множество Мандельброта?


Если коротко: берем точку, на комплексной плоскости. Т.е, точка, это комплексное число, с какими-то мнимой и действительной частью, соответствующей координатам точки на плоскости. Дальше, возводим это число в квадрат, и прибавляем к результату исходную точку. Снова возводим полученное в квадрат, и опять прибавляем исходное число. И так повторяем много-много раз. Смысл вот в чём: исходная точка (число) принадлежит множеству Мандельброта в том случае, если при бесконечно большом количестве итераций вышеприведенного алгоритма, модуль результата не стремиться к бесконечности. Как быстро выясняется, это множество ведет себя весьма и весьма интересным образом, а именно, образует самоподобные (фрактальные) структуры, при приближении к границе множества. Строгой границы множество судя по всему не имеет, поскольку при приближении некоторой точки к этой границе, у нее просто уменьшается скорость стремления к бесконечности, и если предположить, что некоторая точка является границей множества, то всегда найдется другая точка, скорость стремления которой к бесконечности будет меньше, но она, тем не менее, будет обращаться в бесконечность (но с меньшей скоростью, чем предыдущая).


Фрактал Мандельброта в цвете
Рис. №295. Фрактал Мандельброта в цвете

Визуализация множества, по классической схеме, строиться на нахождении количества итераций нашего алгоритма, при котором рисуемая точка обращается в бесконечность. В качестве этой бесконечности берется число 2, ибо доказано, что всякое комплексное число, модуль которого больше 2-х, не принадлежит множеству Мандельброта, т.е. неминуемо обращается в бесконечность. Другими словами, для каждой точки на нашей будущей картинки, мы сначала определяем ее координаты на комплексной плоскости, запоминаем их, потом, в цикле, возводим это комплексное число в квадрат, и прибавляем оригинальные координаты, которые запомнили. Так делаем до тех пор, пока модуль нашего результата будет меньше чем 2. Полученное число итераций можно использовать для задания цвета нашей точки. Необходимо отметить, что количество проводимых итераций нужно ограничивать, поскольку нужно помнить, что если точка принадлежит множеству, то это количество итераций для нее будет равно бесконечности, т.е. мы получим бесконечный цикл. В самом простом случае, ограничение можно сделать равным 255, и просто присваивать его какой-нибудь RGB составляющей выбранной точки, в этом случае мы получим картинку в градации какого-либо цвета. К примеру, в градации серого (все три компоненты цвета равны количеству итераций), картинка классического множества будет выглядеть примерно вот так:


Множество Мандельброта в градациях серого цвета
Рис. №296. Множество Мандельброта в градациях серого цвета

Приведу здесь немного кода на Delphi, с помощью которого была получена эта картинка:


Код: Delphi
type

  //тип комплексного числа, собственно
  TComplex = record
    re,im:Double;
  end;

  //указатель на комплексное число
  PComplex = ^TComplex;

implementation

//базовые операции с комплексными числами

function complex_abs(z:PComplex):Double;
begin
  Result:=sqrt(z.re * z.re + z.im * z.im);
end;


procedure complex_mul(dest, src:PComplex);
var re:Double;
begin
  re := dest.re * src.re - dest.im * src.im;
  dest.im := dest.re * src.im  + dest.im * src.re;
  dest.re := re;
end;

procedure complex_add(dest, src:PComplex);
begin
  dest.re := dest.re + src.re;
  dest.im := dest.im + src.im;
end;

procedure complex_sub(dest, src:PComplex);
begin
  dest.re := dest.re - src.re;
  dest.im := dest.im - src.im;
end;

//возведение в произвольную степень, с помощью тригонометрической формы
procedure complex_powr_trig(dest:PComplex; pow:Double);
var z:TComplex;
zabs,zpow,phi:Double;
begin

  zabs := complex_abs(dest);

  if(zabs = 0)then begin
    exit;
  end;


  zpow := Power(zabs, pow);

  phi := ArcTan(dest.im / dest.re) * pow;

  dest.re:=zpow * cos(phi);
  dest.im:=zpow * sin(phi);
end;


procedure complex_sqr(dest:PComplex);
var re:Double;
begin
  re := dest.re * dest.re - dest.im * dest.im;
  dest.im := dest.re * dest.im  + dest.im * dest.re;
  dest.re := re;
end;



//возвращаем n
function fractal_mandelbrot(c:PComplex):Longint;
var w:Longint;
z:TComplex;
izf:Double;
begin
  z.re:=0;
  z.im:=0;

  w:=254;

  for Result:=0 to w do begin
    //возводим в квадрат
    complex_sqr(@z);

    //прибавляем оригинал
    complex_add(@z, c);
    
    //выясняем, больше ли модуль, чем нужно?
    //число 2000 множно заменить на 2, а можно и не менять
    if(complex_abs(@z) > 2000)then begin
      exit;
    end;

  end;
end;

Остальная часть кода, в простейшем случае, это примерно такой цикл:

Код: Delphi
procedure TForm1.Button4Click(Sender: TObject);
var n,iy,ix:LongInt;
z:TComplex;
begin

  for iy:=0 to (1024 - 1) do begin
    for ix:=0 to (1024 - 1) do begin

      z.re:=(ix - 512) / 256;
      z.im:=(iy - 512) / 256;

      n := fractal_mandelbrot(@z);

      if(n > 254)then begin
        Form1.Image1.Canvas.Pixels[ix,iy]:=0;
      end else begin
        Form1.Image1.Canvas.Pixels[ix,iy]:=rgb(n, n, n);
      end;


    end;
  end;
end;

Числа 1024, 512, и 256 зависят от размера изображения, на вашей форме. В этом алгоритме, выполняется проход по комплексной плоскости от -2 до 2 по осям X и Y, что собственно и приводит к генерации картинки множества Мандельброта. В любом случае, если вы хотите написать программу для более подробного исследования этого множества, то вам придется писать куда более сложный код, поэтому как работает код выше, я думаю, разберетесь уже самостоятельно. В конце статьи, кстати, я приведу полные исходники своей программы.


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


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

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


Далее, несколько слов об оптимизации. Сразу надо сказать, что алгоритм этот очень медленный, и прорисовка областей фрактала в высоком разрешении может занимать от нескольких секунд, до нескольких минут/часов (при высокой глубине). Поэтому, для удобной навигации по множеству, можно сделать “приближенную” прорисовку, т.е. прорисовывать к примеру каждый 16-й пиксель, а в полном разрешении, рисовать картинку только по требованию. Впрочем, не нужно забывать и об оптимизации самого алгоритма, попробую разделить все по пунктам:



1. Количество арифметических операций можно сокращать, например, определение, является ли модуль числа меньше 2:


\[ \sqrt{x^2 + y^2} < 2 \]


можно переписать так:

\[ x^2 + y^2 < 4 \]


и этим, сэкономить на одной операции извлечения квадратного корня.

2. Всё, что разнесено в примере по функциям, можно делать без использования функций, что естественным образом экономит на call и ret командах.


3. Многие функцию фрактала делают рекурсивной. Этого делать не нужно, ибо вариант без рекурсии, с циклом, как выше, работать будет конечно-же быстрее. Ну, конечно же в том случае, если компилятор не обладает способностью к оптимизации хвостовой рекурсии. Но, как я думаю, компилятор имеющийся в Delphi 7 такого делать не умеет. Поэтому, это нужно иметь ввиду.


4. Непосредственно графическая часть: закрашивание пикселей на картинки, естественно, нужно делать не через Canvas.Pixel (что приводит к вызову gdi32.SetPixel() для каждого пикселя), а хотябы с помощью прямого редактирования DIB версии изображения в памяти, т.е. использовать TBitmap.ScanLine[] для получения указателя на конкретную строчку изображения. Как это делать, вы сможете найти в примере, в самом конце.


5. Использовать SSE и/или SSE2 для операций с комплексными числами. При особой изворотливости, вычислять несколько точек одновременно.


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


7. Наконец, в совете, где я говорил о рисовании в 2 прохода, в качестве первого прохода можно брать предыдущий отрисованный кадр, как собственно это и делается в моем примере.


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


Одно из интересных образований фрактала Мандельброта
Рис. №298. Одно из интересных образований фрактала Мандельброта

Теперь попробуем привести основные красоты, найденные мной в нашем множестве. В инфо о изображении я старался указывать координаты, но строго подобных картин в вашей программе скорее всего не получиться получить, по причине к примеру разных алгоритмов вычисления цвета. Одно из интересных образований вы видите выше. Его координаты: Re: -1.175477680975, Im: -0.000000117014, Зум: 4792706. Далее, еще немного о нем:


Без названия
Рис. №300. 2.png

Re: -0.224916283825, Im: -1.116262122278, Зум: 75367985484. Как мы можем заметить, координаты совсем другие. Однако, сама фигура очень похоже на предыдущую. Но не совсем.


Без названия
Рис. №301. 3.png

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


Прибавляя зуму
Рис. №302. Прибавляя зуму

Прибавляя зуму 2
Рис. №303. Прибавляя зуму 2

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


Далее, еще одно интересно образование:


Без названия
Рис. №304. 14.png

Без названия
Рис. №305. 22.png

Без названия
Рис. №307. 18.png

Оно же, только в другой расцветке, и видимо с совершенно другими координатами:


Без названия
Рис. №306. 17.png

Другие наиболее удачные кадры при исследовании множества:


Без названия
Рис. №308. 70.png

Без названия
Рис. №309. 76.png

Биоморф натуральный
Рис. №312. Биоморф натуральный

Это изображение мне очень напомнило живую клетку…


Без названия
Рис. №313. 85.png

Без названия
Рис. №314. 86.png

Без названия
Рис. №315. 88.png

Теперь, приведем несколько инструкций, как найти два основных образования в классическом, не модифицированном множестве, показанных на фото №298 и №306:


Как найти основные фигуры
Рис. №318. Как найти основные фигуры

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


Модифицируем множество Мандельброта.


Но, на самом деле, не только в этом. Красот множества мне показалось мало, и я решил немного поэкспериментировать над основной формулой. Первое что пришло в голову, возведение комплексного числа не во 2, а в произвольную степень. На первый взгляд, это приводило к интересным результатам…


В первые секунды после того, как мысль о произвольных степенях посетила мою голову, я скопипастил функцию возведения комплексного числа в квадрат, потом еще раз… И в итоге, получил примерно это:


Чётные степени
Рис. №321. Чётные степени

Как видим, при повышении степени, в множестве увеличивается количество основных “лучей”, причем это их количество всегда на 1 меньше степени, в которую число возводится, и как опять-же можно заметить, эти лучи всегда располагаются на равных углах друг от друга, тут даже можно провести аналогию с извлечением корня n-ной степени из комплексного числа. Кстати один из способов, которым я возводил комплексные числа в степень, как раз была формула Муавра, требующая представления числа в тригонометрической форме, и как следствие жутко медленная, и вот что интересно: с помощью этой формулы, получалось числа возводить только в чётные степени: 2, 4, 6, 8, …, при попытке возведения числа к примеру в 3-ю степень, получалось примерно это:


Неудачная третья степень
Рис. №323. Неудачная третья степень

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


Неудачная 2.52-я степень
Рис. №326. Неудачная 2.52-я степень

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


Нечётные степени
Рис. №322. Нечётные степени

Здесь степени ведут себя точно также, количество лучей всегда равно числу на единицу меньшему степени… Кстати и надо отметить, что особенно новых форм в степенях отличных от двойки я так и не нашел. Вернее нашел, но… В общем, вот:


Голодная 3-я степень
Рис. №325. Голодная 3-я степень

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


Изменение нулевого Z
Рис. №319. Изменение нулевого Z

Заключается метод в задании Z0, отличном от нуля, и не выходящим за рамки -2 ÷ 2, вернее комплексного числа, с произвольной мнимой и действительной частью, модуль которого не выходил за рамки -2 ÷ 2. Впрочем, метод тоже не ахти какой плодотворный, однако…


Измененный Z0 в третьей степени
Рис. №341. Измененный Z0 в третьей степени

С помощью этого метода, и 3-ей степени комплексного числа, получилась вот такая картинка. Надо таки должное отдать, формы достаточно интересные, но я вам так скажу, оригинальное множество тоже в сторонке не стоит и на месте не дремлет:


Первые шаги к сглаживанию
Рис. №329. Первые шаги к сглаживанию

Классический Мандельброт. Сглаживание 8x
Рис. №330. Классический Мандельброт. Сглаживание 8x

Ещё одна живая клетка
Рис. №328. Ещё одна живая клетка

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


Антиалиасинг


Представленные фото, как раз таки, имеют 8-кратную степень сглаживания, что ощутимо повышает их положительное влияние на человеческий глаз. Кстати, а вот и небольшой, простенький примерчик, наглядно иллюстрирующий работу такого сглаживания:


Антиалиасинг: пример простой
Рис. №331. Антиалиасинг: пример простой

Как видно, даже двукратное сглаживание ощутимо повышает качество картинки, и делает ее более приятной для глаза, но вот разница между 8 и 32-кратным сглаживанием заметна только с большим трудом. Хотя если увеличить прямые линии, на которых хороша видна “лесенка”, то там и эта разница будет ощутимо заметна. Я думаю, что антиалиасинг очень бы ощутимо увеличивал качество таких фракталов как треугольник Серпинского, или кривая Коха, но о них, к сожалению, мы рассуждаем не в этой статье.


А ещё, за всё в этом мире надо платить, и именно поэтому, сложность алгоритма сглаживания неумолимо стремиться к квадрату коэффициента, что приводит к щекотливой ситуации, в которой рендеринг картинки выше с 32x сглаживанием занял ровно минуту, и наглым образом сожрал почти 400 мегабайт оперативной памяти. И это при том, что здесь используется алгоритм практически на чистом SSE2, работающий чуть ли не в 10 раз быстрее алгоритма раннего, которым получены картинки в начале статьи. И это при том, что рендеринг той же картинки без сглаживания вообще, занимал 60 миллисекунд. Но это ещё ладно. Рендеринг трёх картинок выше, всего-лишь с 8-кратным сглаживанием занимал по 10-15 минут, и съедал более 500 мегабайт памяти. Почему так происходит догадаться совсем не сложно, ибо для того, чтобы получить сглаживание по осям X и Y в 2 раза, надо фактически прорисовать 4 исходных картинки, и пикселям картинки-результата, присвоить среднее значения каждых соответствующих 4-х пикселей в картинке высокого разрешения. И получим 2-кратное сглаживание. Нетрудно догадаться, что для 4-кратного сглаживания потребуется прорисовка уже 16 картинок, для 8-кратного 64 картинок, для 16… Считайте, в общем, сами, а то мне как-то, плять, нехорошо даже делается при таких мыслях. Одно только радует, что если сделать обработку фрактала по блокам, а не одной огромной картиной, как реализовано сейчас, то в принципе, проблему пожираемых гигабайтов памяти можно спокойно разрешить. Можно даже сделать распределённое вычисление картинок, на сотне-другой машин, на которых еще будет по паре-тройке процессоров… Стоп, мы замечтались. Просто занимаюсь я фракталами уже 3 дня, и порой кажется, что прозанимаюсь ими ещё минимум год, уж больно интересная тема для размышлений. Ладно. Приведём еще один, более сложный пример влияния сглаживания на картинку, и после вернемся к предыдущим рассуждениям:


Антиалиасинг: пример сложный
Рис. №332. Антиалиасинг: пример сложный

Здесь более наглядно становится заметно влияния сглаживания на качество. А также, на время, необходимое на рендеринг. В частности, пока рендерилась последняя картинка, я успел попить чай…


Впрочем, не всё так плохо, поскольку не стоит забывать так называемые GPGPU технологии, к коим относятся CUDA и OpenCL, которые, кажется, уже упоминались где-то выше. С их помощью, производительность рисования таких веществ может взлететь до небес… Если вдруг мне приспичит заняться этим, то я обязательно напишу об этом статью.


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


Астероиды...
Рис. №333. Астероиды...

Если вам нравятся астероиды, то думаю оцените. И это только один из неудачных экспериментов. Было много других… Например, возведение в степень, зависящую от номера текущей итерации, приводило к чему-то подобному:


Рыбка n-ной степени
Рис. №334. Рыбка n-ной степени

Результатом экспериментов с логарифмом, экспонентой, синусом, и чьей-то матерью, стало что-то подобное:


Адский беспорядок
Рис. №335. Адский беспорядок

Из десятой если не ошибаюсь степени получалось вот что:


Десятая степень
Рис. №336. Десятая степень

Ну и наконец последней каплей для меня стало вот это:


Печально...
Рис. №337. Печально...

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


Ещё один Мандельброт
Рис. №338. Ещё один Мандельброт

Судя по всему, эксперимент
Рис. №339. Судя по всему, эксперимент

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


Множество Мандельброта можно сравнить с бесконечной вселенной, со звёздным небом, особенно это ощущается, когда смотришь на свои координаты в множестве: 12 знаков после запятой, с увеличением в несколько миллиардов раз… Романтика. Других слов и не подобрать. Наверно это и вдохновило написать как эту статью, так и программу… Кстати, ещё несколько слов о программе. Алгоритм генерации множества был мной сначала переписан на SSE2, причем так, что в узком месте цикла расчета каждой точки исключались любые обращения к памяти. Вернее, они вообще исключались практически во всём алгоритме, и все это вместе привело к увеличению производительности примерно в 8 раз. Но мне этого показалось мало. Просидев ещё часик, мной были раскурены некоторые новые инструкции сначала из SSE3, в частности ADDSUBPD, которая взяла на себя половину работы по перемножению комплексных чисел, а затем и DPPD из набора SSE4, явившаяся ни чем иным, как полноценным скалярным произведением, сократившая количество команд, рассчитывающих модуль комплексного числа, ровно до двух штук. Единственное, чего эта команда не сделала, это не налила шампанское, и не заказала девочек. Далее еще были выпилены некоторые лишние команды в циклах, ещё чуть чуть оптимизирован код возведения в квадрат, и все это в совокупэ подняло производительность еще на 10-20%. Но мне и этого показалось мало, и поэтому в ближайших планах сделать возможность рисования в несколько потоков. Но сначала, сделать возможность гибко редактировать цветовую палитру. Думаю, что если это удастся, то можно будет генерировать фракталы такой красоты, коей они нам и не снились… Впрочем, о программе уже точно будет в следующей статье, где, кстати, будут и её исходники. На сегодня всё. Спасибо за внимание! :)

Ключевые слова математика, фрактал, эксперимент, Множество Мандельброта, Delphi 7, комплексное число, бесконечность, SSE, оптимизация, сглаживание, Антиалиасинг
QR-код
Интересное
Как обычно, был семинар по векторной алгебре. Почему-то он показался мне смертельно скучным, может потому что тема плоскостей и прямых в пространстве мне уже была относительна знакома, а может потому что я забыл дома задачник… Впрочем скука продолжалась недолго

Читать »»
Случайные фото