Множество Мандельброта: Рисуем его правильно

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

Впрочем, правильно или не правильно — вопрос скорее религиозный. Речь сегодня пойдёт о той самой программе, которую, некогда, я за несколько дней написал для отрисовки множества Мандельброта. Если вы не читали предыдущую статью, то рекомендую это сделать, дабы узнать, с чего же всё начиналось. В прошлый раз мы рассмотрели простейшую технологию антиалиасинга (то-есть, сглаживания), некоторые приёмы модификации множества, а так-же привели простейший код, позволяющий нарисовать основную фигурку множества. Возможно вы заметили, что код этот не сильно отличается производительностью, особенно если попытаться с его помощью рисовать картинки побольше. Это будет первая проблема, которую мы попробуем сейчас решить…


Проблема первая. Учимся применять SSE.


Оптимизация, это вообще такая вещь, про которую можно говорить сколь угодно долго, придумывать всё новые и новые ухищрения дабы сэкономить еще парочку наносекунд… Впрочем то, что удалось реализовать мне, производительность прорисовки множества поднимает почти в 6 раз, а производительность накладывания на него цветовой палитры — примерно в 8-10 раз, по крайней мере после такой оптимизации стало возможно наблюдать результат в реальном времени, двигая мышкой ползунки цвета в программе. Попробуем привести конкретные скриншоты, показывающие время отрисовки:


Результат оптимизации, или торжество SSE
Рис. №883. Результат оптимизации, или торжество SSE

Надо сказать, что шестикратный прирост производительности — результат, очень мягко говоря, неплохой. Как вы уже заметили, достигнут он с помощью расширения SSE, причём используются в нём команды сразу из SSE, SSE3 и SSE4. Фундаментальный минус такого подхода — в нём очень сложно как-либо изменить алгоритм генерации множества, во первых потому что для каждой, вроде-бы простой последовательности вычислений, приходится придумывать эквивалент на SSE, который далеко не всегда бывает столь-же очевиден по функциональности (в большинстве случаев наоборот, код с использованием SSE совершенно не читабелен), а во вторых — при таком подходе мы имеем, грубо говоря, всего 8 переменных — именно столько в x86-й архитектуре XMM регистров. Здесь я иногда завидую обладателям x86-64 — у них этих регистров уже целых 16. Нет, конечно-же, количество переменных можно “увеличить”, если в качестве переменных использовать память, но я думаю очевидно, что за обращение к памяти в подобных алгоритмах приходиться очень и очень дорого платить…


Давайте мы немного расскажем, что-же такое SSE и почему я его так люблю. Если коротко, то это небольшое расширение современных процессоров архитектур x86 и x86-64, реализующее принцип SIMD, т.е. “Одна Инструкция, Много Данных”, и конкретно, позволяющее производить арифметические (и не только) действия одной командой, над несколькими числами (операндами) одновременно. Имеет несколько версий, и даже в некоторой из них расширяет функциональность MMX (подобное SSE расширение, только старее, и предназначенное для работы с целыми числами). Расширение SSE использует 8 (или 16 в x86-64) специальных, 128-битных регистров, именуемых XMM-регистрами, позволяющих, как быстро можно заметить, производить действия одновременно над 4-мя числами типа float или long int (4 байта, SSE1 и расширения MMX), 2-мя числами double (8 байт, в основном SSE2), и сразу 8-мью однобайтовыми величинами (некоторые команды MMX)


Список команд и возможностей SSE весьма впечатляет. Можно складывать, вычитать, умножать, делить, извлекать корни, делать сдвиги, сравнивать, и всё это над несколькими парами операндов одновременно. На лицо довольно основательный потенциальный прирост в скорости… Ещё в SSE есть интересные команды для быстрой работы с памятью, например MOVNTPS. Кстати работа с памятью в SSE производиться только по выровненным (кратным 16) адресам, по причине значительного улучшения быстродействия по сравнению обычными адресами. Это порождает как плюсы, так и геморрой… Есть конечно команды, работающие с не выровненными адресами (например, MOVUPS), но использовать их не рекомендуется — ибо медленно. DRAM память так устроена, что адреса кратные 16 как раз показывают на начала специальных ячеек, которые считываются за 1 такт, и поэтому не выровненный адрес всяко заденет 2 ячейки, и такта придётся ждать уже 2… Впрочем всё это выдержки из матчасти. Подробнее про SSE можно почитать, например, в Википедии.


Вернёмся теперь к нашему барану, и рассмотрим, как же мы из 5 с половиной секунд получили чуть меньше одной секунды. Для этого приведём сначала сокращённый вариант старого кода, который, собственно, использовался в не оптимизированной версии программы. Мы его уже приводили в предыдущей статье, но приведём еще раз:


Код: Delphi
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);
    
    if(complex_abs(@z) > 8)then begin
      exit;
    end;

  end;
end;

procedure Draw();
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;

Можно сразу заметить, что сложность сего алгоритма составляет примерено 
\[ O(H \cdot W \cdot D) \]
, то-есть зависит от произведения высоты, ширины, и глубины прохода, что само по себе, довольно печально. Но с этим мы, похоже, поделать ничего не сможем. Хотя если подумать… Мы только можем сказать, что функция fractal_mandelbrot() в данном коде является самым узким, критичным по времени местом, вернее не она сама, а цикл в ней. Первая идея — от функции избавиться в принципе, то-есть код её вставить в то место, откуда она вызывается (заинлайнить), чем можно сэкономить примерно W ⋅ H мусорных команд call и ret. Сэкономим пару миллисекунд в лучшем случае. Дальше, можно избавиться от лишней операции взятия квадратного корня при вычислении модуля комплексного числа, как я уже где-то писал. Просто взгляните на тождество:

\[ \sqrt{x^2 + y^2} < 2 \Leftrightarrow x^2 + y^2 < 4 \]
Капитан Очевидность заявляет

Дальше, уже не всё так просто. Ну, естественно, можно избавиться от функций complex_*, которыми пестрит код выше, заменив их, эквивалентным кодом. Опять сэкономим миллисекунды, правда на этот раз экономия будет уже порядка W ⋅ H ⋅ D, что довольно значительно. Но дальше… А вот дальше, нам придётся опускаться на уровень ниже, и начинать мыслить уже ассемблерно. И вот здесь, нам на помощь, уже приходит SSE. Но, не сразу…


Отладчик Delphi 7 и MMX
Рис. №884. К слову, Delphiйский отладчик не умеет показывать XMM регистры...

Для начала, задумаемся, на одну маленькую тему, а именно, где у нас, собственно, хранятся все переменные, используемые программой. Это, конечно же очевидно, что в стековой памяти, в бо́льшем своём числе. Ключевое слово здесь — в памяти… Обращение к которой, в таких критических местах как наш третий цикл — операция, по сути говоря страшно расточительная. Если посчитать, сколько раз происходит обращение к памяти на запись/чтение, то даже тот довод, что в цикле этом по любому будет использоваться кэш процессора, уже не катит. Третья наша основная цель — ПОЛНОСТЬЮ исключить обращение к оперативной памяти при проведении вычислений. Серьёзно, вообще полностью. Все данные хранить будем только в регистрах процессора. Где-то я читал, как один автор реализовал этот алгоритм, храня данные в стеке математического сопроцессора (стеке FPU/MMX), и это безусловно было круто, но уже, к сожалению, не модно. Мы пойдем дальше, и данные будем хранить… Да! Вы несомненно угадали! В XMM регистрах! И для подавляющего числа операций мы будем использовать возможности, собственно говоря, SSE.


Итак, первым делом определимся с точность. Комплексные числа мы будем хранить как и ранее, в виде двух чисел double, поскольку при вычислении фрактала с очень большим увеличением, точность начинает играть далеко не последнюю роль. А ещё, одно комплексное число будет удобно размещаться в одном 16-ти байтовом регистре, что, несомненно, довольно удобно. В нашем распоряжении есть 8 регистров XMM, 7 регистров общего назначения (4-х байтовых, исключая ESP), и в запасе ещё стек FPU, с его 8-ю регистрами, которые по совместительству, кстати, ещё являются регистрами MMX. Последние, правда, оказались совершенно не нужны.


Дальше, нужно разобраться с константами. Можно насчитать 5 констант: 1024, это разрешение единичного квадратика комплексной плоскости при масштабе 1:1, 4 (было заменена на 8) — с ней будет сравниваться модуль комплексного числа, текущее положение на комплексной, коэффициент масштаба, который будет означать масштаб… Хотя в принципе, в него можно включить и константу 1024, поэтому её уже не считаем. Ещё нужно будет хранить следующие вещи:

1. Указатель на область памяти, в которую мы будем складывать готовые вычисленные значения;
2. Текущее значение комплексной переменной C при вычислении, собственно, фрактала;
3. Чуть позже выясняется — где-то надо хранить 4 вычисленных значения, прежде чем отправить их в память;
4. Собственно, счётчики цикла ширины и высоты, для них сойдут ECX и EDX (в итоге, остался только ECX).

Ну… И пожалуй всё. Остальное придёт уже в самом процессе. Думаю теперь мы готовы увидеть, что-же в итоге получилось… Нужно сразу предупредить, код весьма сложный. Начнём с начала:


Код: Delphi
//ускоряемся, ибо уже как-бы пора
//функция может работать только с картами, с шириной кратной 4-м
//выравненную память мы готовим в другом месте
procedure fractal_mandelbrot_fast(map:PTraceMap; nmax, w,h, wu,hu:LongInt;
  xstart,ystart, zscale:Double) stdcall;
var
cconst:Double;
t0, t1:Double;
begin
  cconst:=8.0;

  Form1.ProgressBar1.Position:=0;
  Form1.ProgressBar1.Max:=h;

  t0 := perf_time();

  asm
    {$I ./inc/mandelbrot_fast.asm}
  end;

  t1 := (perf_time() - t0) * 1000;

  Form1.LabelTime.Caption:=Format('%0.2f ms', [t1]);
  Form1.ProgressBar1.Position:=0;
end;

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


Код: Delphi
procedure aligned_get_mem(var p1,p2:Pointer; count, edge:LongInt);
begin
  GetMem(p1, count + edge);
  p2:=Pointer((DWORD(p1) + edge) and (not (edge - 1)));
end;

Выглядит таинственно и непонятно, но мы поясним. Функция выделяет немного больше памяти, чем нужно, а именно больше на размер выравнивания, которого нужно достичь, в один указатель записывает оригинальный адрес на выделенную память, а в другой — ближайший адрес, кратный 16-ти. Последний получить довольно просто, прибавив к оригинальному адресу 16, и сделав логическое AND с NOT (16 - 1) (подумайте, почему так?). 16 здесь — требуемая граница выравнивания, и есть ограничение: граница эта должна быть степенью двойки, то-есть состоять только из одного бита. В целом, решение это показало себя отлично.


Ещё одна непонятная функция, perf_time(), занимается ни чем иным, как подсчётом прошедшего времени с помощью специальных Windows API:


Код: Delphi
function pref_snap():Int64;
begin
  QueryPerformanceCounter(Result);
end;

function perf_time():Double;
begin
  Result:=pref_snap() / pref_freq;
end;

В pref_freq здесь должно быть до этого записано значение, возвращённое функцией QueryPerformanceFrequency(). Если кто-то не знает матчасть, то функция QueryPerformanceCounter() возвращает количество тиков абстрактного счётчика с момента старта системы с некоторого опорного значения, (я сам видел, как венда переписывала при загрузке регистр MSR_TSC), а QueryPerformanceFrequency() возвращает количество таких тактов, происходящих за одну секунду. Соответственно, чтобы узнать нам, сколько секунд прошло с того опорного момента, нужно текущее значение счётчика поделить на его частоту. Ну, а когда мы научились с высокой точностью фиксировать определённые моменты времени, нам совершенно не составит труда замерять скорость работы некоторых кусков кода, как собственно это и делается в коде функции fractal_mandelbrot_fast() выше. Надо сказать ещё, что точность таких замеров весьма высокая, и позволяет фиксировать промежутки времени с точностью до долей микросекунд. Аппаратная реализация абстрактного счётчика может быть разной, где-то этим занимается микросхемка времени на материнской плате, а например моя система использует для замеров команду RDTSC, которая возвращает количество тактов, совершенных процессором с момента старта, или последнего обнуления регистра MSR_TSC. Учитывая, что например у меня частота CPU составляет 3.2 GHz, нетрудно заметить, что в теории, с помощью RDTSC можно отмерять промежутки времени в пределах долей наносекунд… Правда на практике этим заниматься особенно не следует, ибо процессор часто на одних и тех-же командах ведёт себя по разному. Суперскалярность же…


Кстати, суперскалярность, это ещё одна такая штука, которую нам бы неплохо было бы использовать в нашем алгоритме. Надо сказать, что сделать это почти не получилось… Суперскалярность, это вообще когда процессор начинает некоторые команды в вашем коде выполнять почти параллельно, в том случае, если они как-бы не зависят друг от друга, то-есть используют разные регистры, разные области памяти и прочее. Позволяет значительно (у меня получалось с 2.5 мкс довести код до 2.0 мкс) поднять производительность некоторых программ. Наш код, как выяснилось, почти не содержит областей, которые можно таким образом распараллелить… То, что получилось в итоге, и без того было весьма сложно в понимании и реализации…


Ладно, оставим философию. Пришло время сорвать покровы, и наконец-то увидеть наш код, вернее, листинг подключаемого файла mandelbrot_fast.asm:


Код: Delphi
//алгоритм должен быть ЭКСТРЕМАЛЬНО быстрый
//поэтому, он будет использовать вперемешку SSE и SSE2 и максимальную оптимизацию
//надеюсь, что получится алгоритм сделать суперскалярным (к сожалению, почти не получилось)
//в том числе ассемблерные циклы, и movntps
//стек будет использоваться всего лишь h раз
//не CUDA конечено, но... прирост скорости в 6 раз! эгегееее!!

//сохраняем регистры, ибо мы их сейчас будем портить
push edi
push esi
push ebx
push ebp

//инецылезируем константы
mov eax, w
imul eax, h
shl eax, 2

//указатель карты должен показывать на последнее число
//ибо для того чтобы выходным резудбтатом можно было пользоваться, мы его будем на каждой итерации уменьшать
//поскоку основная проблема сдешних циклов, это декремент ecx с проверкой на ноль, что не совсем соотвествует
//нормальному человечьему мышлению, когда счетчкики, как у нормальных людей, увеличивваются, а не наоборот..
//и получается, что прорисовку алгоритм делает С НИЗУ, и а не сверху
//ну блять, ну вот все через жопу, все не как у людей! :))
mov edi, map
add edi, eax
sub edi, 16

//запоминаем переменные
mov ebx, nmax
mov eax, w

//во все-х случаях, в младших разрядах будет храниться x, а в старших y
//надо этого всегода придерживаться и учитывать, ибо movaps [структура], регистр вернет именно такой порядок..

//вычисляем коэффициенты масштабирования, а также начальные значения
//для всего этого изобилия, у нас будут использоваться регистры xmm6 (коэффициенты) и xmm7 (смещения)
movlpd xmm6, ystart
shufpd xmm6, xmm6, 1
movlpd xmm6, xstart


//грузим коэффициенты размерности...
cvtsi2sd xmm7, hu
shufpd xmm7, xmm7, 1
cvtsi2sd xmm7, wu

//надо всю эту хрень умножить(вернее поделить) на масштаб
//грузим масштааааб
movlpd xmm0, zscale
shufpd xmm0, xmm0, 0

//ии деелииимм! (вернее множим)
mulpd xmm7, xmm0
//и вот теперь мы имеем нужные коэффициенты)

//еще один момент: в xmm5 пихаем константу, отвечающую за максимальный модуль
//xorpd xmm5, xmm5
//movlpd xmm5, ss1
//shufpd xmm5, xmm5, 1
movlpd xmm5, cconst

//наконец, xmm4 будет содержать промежуточное значение итераций, 4 штуки
//сдвиг будем выполнять с помощью shufps
//и не вздумайте ляпнуть мне тут, что жто невозможно!!
xorpd xmm4, xmm4
//просто кэп вам наменул, есть какието более человечные команды логического сдвига, из MMX..
//товарищь Кэп, мы как нибуть без вас разберемся, ок?!1

mov ecx, h

//первый цикл по строчкам
@for_height:

  //запоминаем счетчик
  mov esi, ecx
  mov ecx, eax

  //и закручиваем второй цикл, уже по пикселям
  @for_width:
	//тестирование комплексного числа должно производиться только в памяти, и нигде более
	mov ebp, ecx  //сохраняем счетчик цикла
	mov ecx, ebx  //устанавливаем nmax

	//далее нужно вычислить координаты точки
	//вычисляем... по аналогии с тем что выше
	//только мы помним, что нужные значения находятся в регистрах!
	cvtsi2sd xmm0, esi
	shufpd xmm0, xmm0, 1
	cvtsi2sd xmm0, ebp

	//множим на необходимые коэффициенты
	divpd xmm0, xmm7

	//прибавляем стартовое положение
	addpd xmm0, xmm6

	//вроде все готово. xmm0 содержит константу c. xmm6 и xmm7 заняты, помним!
	//и первым его элементом будет xmm1
	//в котором, кстати, изначально ноль
	xorpd xmm1, xmm1


	@trace_step:

	  //возводим эту хрень в квадрат

	  // ; re := dest.re * dest.re - dest.im * dest.im;
	  // ; im := dest.re * dest.im + dest.im * dest.re;

	  movapd xmm3, xmm1

	  //у второго регистра меняем все местами
	  //movapd xmm2, xmm1
	  shufpd xmm3, xmm3, 1

	  //перемножаем (действие первое)
	  mulpd xmm3, xmm1  //; re*im + im*re
	  mulpd xmm1, xmm1  //; re*re - im*im

	  //xmm1 : re*re - im*im
	  //xmm3 : re*im + im*re


	  //теперь, нужно у первого регистра старшую часть отнять у младшей
	  //а у второго обе части сложить
	  //и сделать так, чтобы результаты стали частями третьего регистра
	  //и тогдааа!!! мы получим то что хочим, тоесть наше комплексное прозиведение

	  movapd xmm2, xmm1

	  unpcklpd xmm1, xmm3
	  unpckhpd xmm2, xmm3

	  //в SSE3, оказывается! еще есть команда ADDSUBPD
	  //и очень похоже, что она - именно то что нам нада
	  //ADDSUBPD XMM1, XMM2
	  db $66, $0f, $d0, $ca
	  //ИЧСХ, Работает, Сука!


		  //хорошая попытка сделать быстрее. к сожалению, неудачная
		  (****************************************
		  //И... Снова SSE3! :)
		  //грузим числа в хитром виде
		  //movddup xmm3, xmm1
		  db $f2, $0f, $12, $d9
		  shufpd xmm1, xmm1, 1
	
		  //movddup xmm2, xmm1
		  db $f2, $0f, $12, $d1
	
		  //множим
		  mulpd xmm3, xmm1
		  mulpd xmm2, xmm1
	
		  //меняем числа местами в первой части
		  shufpd xmm3, xmm3, 1
		  ***************************************)

	  //УсТаРеЛо

		  (********************************************************************
		  //теперь в xmm1 младншие части бывших xmm1 и xmm3, а в xmm3 - старшие
		  //остается только их както хитро сложитьвычесть
	
		  xorpd xmm3, xmm3
		  unpckhpd xmm3, xmm2
	
		  //стираем младшую часть и прибавляем
		  //а после этого сохраненную младшую часть вычииитаем!
	
		  addpd xmm1, xmm3
		  subsd xmm1, xmm2
		  //*********************************************************************)

		  //это для дебага
		  //////////////////////////
		  {pushad
		  movaps [edi], xmm1
		  mov eax, edi
		  call memview
		  popad  }
		  //////////////////////////

	  //ну все! теперь в xmm1, кажись, валяетсо комплексный квадрат
	  //а поскоку он там лежит, к нему надо приебtybnm xmm0!
	  addpd xmm1, xmm0

	  movapd xmm3, xmm1

	  //и еще одна есть инструкция, интересная, но уже в SSE4
	  //и умеет она делать скалярное произведение, и зовется она DPPD...
	  
	  //DPPD XMM3, XMM1, 49 	 //1 | 16 | 32
	  db $66, $0f, $3a, $41, $d9, 49
	  //РАБООТААЕЕЕТТТ!!!! мой CPU поддерживает SSE4! АААА!!!
	  //а также, новостью является то, что компейлятор, встроеный в Dev-C++, тоже поддерживает SSE4.
	  //чтож, похвально. чего вот не скажешь о нешей дельфи 7.. 


	  //дальше надо найти модуль по оптимизированной технологии
	  //оптимизированная технология была безнадежно устарела, и заменена на sse4.dppd

		  (********************************************************************
		  movapd xmm2, xmm1
		  mulpd xmm2, xmm2
		  //перемножили, теперь части надо сложить
		  movapd xmm3, xmm2
		  shufpd xmm3, xmm2, 1
		  //переставили, складываем
		  addsd xmm3, xmm2
		  //*********************************************************************)

	  //младшая часть xmm3 теперь содержит квадратный модуль
	  ucomisd xmm3, xmm5
	  //сравниваем, и нас интересует только значение флага CF

	  //если меньше, то мы кажется живы
	  jnc @are_dead

	//как выяснилось, такой цикл работает.. быстрее?!
	dec ecx        
	//test ecx, ecx
	jnz @trace_step

	//loop @trace_step

	@are_dead:

	//тут у нас яяяявноо чему то равен ecx, и чем, нам надо во первых выяснить, 
	//во вторых преобразовать в нормальное, человеческое представление
	mov edx, ebx
	sub edx, ecx
	//ураа!! edx содержит искомео количество итераций! ураа!!
	//признаться, я слегка запарился его искать

	//поскольку известно, мы складываем его в xmm4
	//а сначала его двигаем, блять! блять, блять, и блять еше раз!
	//вроде так...
	shufps xmm4, xmm4, 147    //10,01,00,11
	cvtsi2ss xmm4, edx
	//правда еще есть команды логических сдвигов из MMX..

	//лулз получается в том, что в выходном массиве у нас получаются ФЛОАТЫ, БЛЯТЬ
	//ладно хуй с вами, флоаты так флоаты... учитывая что у нас будет еще линейная 
	//интерполяция, это наверно даже к лучшему
	//просто пока это писалось, никто даже и не догадывался
	//о существаоании интересной такой команды CVTTPS2DQ...

	//смысл дальнейшего кода в сбросе регистра xmm4 по требованию, т.е по достижению его заполненности
	//далее олгаритм довольно прост. во первых, интересуемся, чему равен j:
	//регистр ecx тут уже не нужен, он переписывается несколькими командами ниже
	mov ecx, ebp
	dec ecx
	and ecx, 3 //если кратно 4м
	jnz @skip_flush

	  //сбрасываем этот ебучий регистр
	  movntps [edi], xmm4

	  //двигаем адрес на карте (отнимаем)
	  sub edi, 16

	@skip_flush:
	//весь смысл хрень выше, как вы могли уже догадаться, это сброс 16 байтов в память каждые 4 пикселя
	//на самом деле, работает как часы

	//вспоминаем значение счетчика
	mov ecx, ebp

  //цикл по элементам строки
  dec ecx
  //test ecx, ecx
  jnz @for_width
  //loop @for_width

//цикл по строкам, вспоминаем его щщоттчик
mov ecx, esi

//тут запилим небольшой прогресс, срабатывающий раз в дохрена строк
mov edx, ecx
and edx, $1F

//небольшой костыльчик для прогрессбара
jnz @no_progress
	 
  pushad

  sub esp, 128

  movups [esp], xmm0
  movups [esp + $10], xmm1
  movups [esp + $20], xmm2
  movups [esp + $30], xmm3

  movups [esp + $40], xmm4
  movups [esp + $50], xmm5
  movups [esp + $60], xmm6
  movups [esp + $70], xmm7

	push ecx
	call cb_progress

  movups xmm0, [esp]
  movups xmm1, [esp + $10]
  movups xmm2, [esp + $20]
  movups xmm3, [esp + $30]

  movups xmm4, [esp + $40]
  movups xmm5, [esp + $50]
  movups xmm6, [esp + $60]
  movups xmm7, [esp + $70]

  add esp, 128

  //досрочный выход, если там где-то была нажата ESC
  test eax, eax
  popad
  jnz @lets_fuck

@no_progress:

dec ecx
//кэп вам снова говорит, что test делать не надо, ибо dec УЖЕ выставляет все нужные флаги..
//test ecx, ecx
jnz @for_height

//вспоминаем старое значение регистров
@lets_fuck:

pop ebp
pop ebx
pop esi
pop edi

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


Итак, в качестве счётчика циклов у нас всегда используется регистр ECX. Изначально это происходило потому, что для организации циклов использовалась команда LOOP, но потом выяснилось, что связка команд DEC ECX и JNZ ВНЕЗАПНО работает быстрее, что для меня было жёстким когнитивным диссонансом. Причём выигрыш был значительный, вплоть до десятков миллисекунд. Поскольку циклов у нас тут 3, и во всех используется один и тот-же регистр, его значение явно надо где-то хранить, и делаем мы это нет, не с помощью стека! Для этого у нас, собственно, эксплуатируются регистры EBP и ESI, их значения, кстати, в коде где-то неоднократно используются. Кстати, для того чтобы подружить XMM регистры и обычные регистры x86, используются с успехом команды CVTSI2SD и CVTSI2SS, которые конвертируют 4-х байтовое целое в число с плавающей точкой типа double и float, соответственно, и записывают их в младшую часть указанного XMM регистра. Готично, однако.


Далее, попробуем осилить возведение комплексных чисел в квадрат с использованием SSE. Формула умножения комплексных чисел, вообще, имеет примерно такой вид:


\[ (a+bi)\cdot(c+di)=ac+bci+adi+bdi^2=(ac-bd)+(bc+ad)i \]
Умножение комплексных чисел

Но не сто́ит сразу унывать… Нам нужно вообще-то умножить число само на себя, чтобы собственно возвести его в квадрат. Ну, а формула возведения комплексного числа в квадрат вырождается приблизительно в следующее:


\[ (a \cdot a - b \cdot b) + i(b \cdot a + a \cdot b) = (a^2 - b^2) + i2ab = \underline{(a^2 - b^2) + i(ab + ab)} \]
Возводим комплексные числа в квадрат

Смотреть нужно на подчёркнутую часть формулы, в ней счастье! Просто очевидно, что половину работы можно возложить на 2 команды попарного умножения MULPD, то-есть просто умножить число само на себя, и потом снова умножить само на себя, только перевёрнутое. Перевернуть можно с помощью SHUFPD, например. Этот замечательный алгоритмец и реализуется в вышеизложенном коде… Полное описание команд мы не будем приводить, лучше приведём ссылку на ОЧЕНЬ хороший мануал от самого Intel, в котором всё сверхдосканально разжёвано. Мануал сей, собственно, состоит из нескольких частей, скачать которые можно, например тут (часть первая), и здесь (часть вторая).


Но умножение, это только половина работы. Как видно из формулы, после 4-х умножений, нужно как-то хитро половину результатов вычесть, а половину сложить. Ну, вы поняли. Надо сказать, это одно из мест где моск пришлось немного поломать. Было всего 3 варианта реализации (все есть в коде под комментариями), первый брутальный на SSE2, второй неудачный с помощью команды MOVDDUP, и наконец третий, ставший финальным, использующий специально для этого предназначенную команду ADDSUBPD и связку из UNPCKLPD и UNPCKHPD, собственно говоря как он работает, отражено в коде. Такой вариант оказался самым быстрым, даже выигрывающим у немножко другой комбинации, который я нашел в гугле. Отметить только надо, что команды MOVDDUP и ADDSUBPD появились только в SSE3, и поэтому моя любимая Delphi 7 о их существовании даже не подозревала, из-за чего в коде пришлось писать их в виде бездушных циферок, любезно сгенерированных попавшимся под руку Dev-C++. Кстати то, что Dev-C++ осведомлена даже о самых новых командах из SSE4, для меня было приятной новостью. Ну да ладно…


Следующим MindFuck'ом стал вроде-бы нехитрый алгоритм вычисления модуля нашего комплексного числа, являющегося по совместимости ничем иным, как скалярным произведением вектора, составленного из действительной и мнимой части. Вернее, не совсем модуля — квадратный корень вычислять, как мы помним, совсем не обязательно. Сначала конечно были мною использованы возможности одного лишь SSE2, но когда я почитал о возможностях SSE4, ситуация сразу феерично прояснилась — оказалось, что среди набора его инструкций есть интересная такая парочка команд, DPPS и, подходящая под наш 8-ми байтовый случай DPPD, предназначенная для вычисления ни чего иного как… Скалярного произведения. Да блядь! Одной командой! Ещё и результатом вычисления можно было засрать целый регистр с помощью маски. Естественно, по производительности с этой инструкцией практически любому коду тягаться было бессмысленно. Брутальности ей добавлял угрожающе длинный опкод, занимающий целых 6 байт, в коде выше её поэтому трудно не заметить. Поскольку, кроме скалярного произведения, нам в общем-то ничего нужно и не было, осталось только сравнить посчитанное число с константой 8.0 с помощью UCOMISD, проверить флаг CF с помощью JNC, и…


И наступает, наконец-то, завершающая часть расчёта, а именно сброс посчитанных результатов в память. Мне пришла хорошая идея делать это только на итерациях, кратных 4-м, то-есть сбрасывать сразу по 16 байт, и как можно сразу догадаться, делать это с использованием выровненного адреса, и инструкции MOVNTPS. Правда это повлекло за собой ограничение, по которому ширина строк должна быть кратна 4-м, но это в целом-то и не страшно, ведь память-то у нас всё равно передаётся выровненная… Ну, как определять кратность номера итерации 4-м объяснять я думаю не нужно, а вот сдвиг регистра-накопителя… Да-да, данные ведь, перед тем как отправить в память, надо где-то накапливать, вот для этого у нас и используется XMM регистр, который после каждой итерации циклически сдвигается на 4 байта, подставляя таким образом место для нового float числа, и при полном заполнении образует именно такой симпатичный кусочек данных, который и должен сохраняться в память. Вот так всё хитро :) Заметить надо, что сдвиг реализован с помощью известной уже команды SHUFPS, мощь которой иногда поражает. Вся фишка в используемой маске 147, которая в двоичном представлении выглядит как 10`01`00`11, где хорошо видно, что пары битов образуют сдвинутую убывающую последовательность, а поскольку для команды эти пары являются индексами четвёрок байтов в регистре, так и получается, что регистр просто циклически сдвигается на 1 4-х байтовый элемент. И никакой магии тут на самом деле нету…


Последней доработкой тут была реализация обновления прогресс-бара по мере прорисовки карты, чисто, надо сказать из эстетических соображений. После нескольких трудноуловимых ошибок оказалось, что неплохо бы перед вызовом функции прогресса сохранять все XMM регистры, поскольку иногда, в зависимости от чётности фазы луны, после вызова прогресса, содержимое наших регистров местами сильно портилось. Трудно, конечно, представить, что Windows API или Дельфийский рантайм где-то использует XMM регистры или вообще SSE, но факт, сука, есть факт. Кстати ещё один момент, красноречиво говорящий нам о том, что использовать новые команды в поздних версиях SSE весьма полезно (речь идёт в данном случае о DPPD и ADDSUBPD) — поглядим на говорящие за себя скриншоты:


Сравнение производительности версий SSE
Рис. №888. Сравнение производительности версий SSE

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

1. Полученную карту высот нужно сжать в N раз, где N — глубина сглаживания;
2. На карту нужно наложить цвет, с помощью некоторой функции от высоты в точках.

Сразу надо сказать, сглаживание лучше реализовать максимально просто, а максимально просто, это значит прорисовывать карту в размерах W ⋅ N и H ⋅ N, где W и H — ширина и высота, а N — коэффициент, или глубина сглаживания, затем, для точек результирующей карты, брать глубину, равную среднему арифметическому глубин в квадратах размером N ⋅ N. Я думаю против того, что это максимально просто, никто возражать не будет, и поэтому, сразу приведём очередной огрызок кода, который, собственно говоря, это всё реализует:


Код: Delphi
//размер должен быть строго кратен четырём
procedure fractal_map_generate(map:PTraceMap; mmw,mmh, wres:LongInt; 
  xstart,ystart, zscale:Double; approx:Double; zdepth:Longint);
var
  farx, zmw, zmh, zmz, darx:Longint;
  tmpp, tmpo:PTraceMap;
  i, j, li, lj, jss, jsm: Longint;
  acc: Single;

begin

  if(approx > 1)then begin
    exit;
  end;

  farx := round(1 / approx);

  darx := farx * farx; 

  zmw := Round(mmw) * farx;
  zmh := Round(mmh) * farx;
  //размеры теперь 100% выровнены

  zmz := Round(wres) * farx;

  //отводим память для большой картинки
  //указатель на карту, как не трудно догадаться, показывает на маленькую карту
  aligned_get_mem(Pointer((@tmpo)^), Pointer((@tmpp)^), zmw * zmh * 4, 16);

  fractal_mandelbrot_fast(tmpp, zdepth, zmw, zmh, zmz, zmz, xstart, ystart, zscale);

  //интерполируем
  for j:=0 to mmh - 1 do begin

    jsm := j * mmw; 

    for i:=0 to mmw - 1 do begin

      acc:=0;

      //ищем средние значения квадрата
      for lj := j * farx to j * farx + farx - 1 do begin

        jss := lj * zmw;

        for li := i * farx to i * farx + farx - 1 do begin
          acc := acc + tmpp[jss + li];
        end;
      end;

      map[jsm + i] := acc / darx;
    end;
  end;

  FreeMem(tmpo);
end;

Тут всё просто. Параметр map показывает на карту требуемого размера, внутри функции создаётся новая карта, размерами W ⋅ N × H ⋅ N, в нее генерируется соответствующих размеров фрактал, и после генерации, реализуется, собственно говоря, вышеописанный алгоритм. Конечно, с билинейной интерполяцией он имеет не слишком много общего, но работает более чем прекрасно. В силу небольшой квадратной сложности я решил не использовать SSE при его реализации, а так и оставить: даже для очень больших карт при сглаживании 16× время работы алгоритма редко бывало больше секунды, по крайней мере на моей двухъядерной системе с частотой 3.2 GHz. Кстати двухъядерность тут роли почти не играет, используется для рисования, естественно, только одно ядро. Как сделать так, чтобы можно было распараллеливать процесс вычислений, изложим мы, скорее всего, в следующей статье. А пока что…


Множество Мандельброта...
Рис. №886. Да будет цвет!

Проблема вторая. Да будет цвет!


А пока что, у нас остаётся проблема цвета, которая на самом деле далеко не такая простая, как кажется проблема сглаживания. Вообще, алгоритмы для задания функции цвета от “глубины” можно выдумывать разные, на первых порах такой функцией вообще служил массив, заполненный случайными числами, но по понятным эстетическим причинам, такой школьный подход сейчас нам малоинтересен. Я, впрочем, решил особо далеко не залезать, и в функцию эту воплотить возможность создания градиентов цвета, задаваемых опорными точками, почти как это мы делаем в фотошопчике. Местоположением точки, в данном случае, будет являться значение глубины, ну а значения между точками будут генерироваться с помощью какой-нибудь интерполяции. Первая версия программы реализовывала только линейную интерполяцию, но потом, немного подумав, я решил добавить туда более гладкую косинусную интерполяцию, и как выяснилось — не зря. По крайней мере, на мой субъективный взгляд, более контрастные и сочные картинки, получающиеся с использованием косинуса смотрятся лучше, чем более мутные и блёклые, которые генерирует линейная интерполяция. Поскольку, это уже вопрос религиозных и других пристрастий, в программе метод интерполяции можно переключать. Вот, кстати, неплохая попытка показать разницу двух методов:


Сравнение методов интерполяции
Рис. №887. Сравнение методов интерполяции

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


Рассматривать процедуру генерации таблицы цветов более детально, я думаю, не сто́ит. Просто берём, идём по всем имеющимся точкам, и промежутки между ними заполняем интерполированным цветом. У каждой точки 2 параметра: цвет и глубина, на которой точка находиться, в моей программе координаты задаются последовательно, то-есть все последующие точки зависят от предыдущей. Может кому-то это может показаться не слишком удобным, но на мой субъективный взгляд, метод оказался довольно удачным, учитывая весьма причудливую природу фрактальных карт, и большую непредсказуемость общей картины даже при незначительном смещении всех точек цвета.


Дальше следует пояснить алгоритм выборки цветов из таблицы. Мы здесь натужено вспоминаем мозгом головы, что в карте числа у нас, мягко говоря, не совсем целые, в отличие, например, от индексов в таблице цветов. Как быть? Ответ один — интерполировать. И поскольку здесь довольно-таки важна высокая производительность, мы обойдёмся без косинусов и прочих сложностей, и снова применим линейную интерполяцию. Мы просто будем интерполировать значения двух соседних элементов таблицы цветов, между которыми будет заключён наш не совсем целый пиксель. А ещё, вся соль в том, что делать мы это снова будем на SSE! По той простой причине, что иным способом не получиться в реальном времени перерисовывать весь фрактал, скажем двигая ползунок цвета в какой-то точке. Ведь сложность наложения цвета квадратична, что уже начинает быть слегка критично по времени, в отличие, например, от алгоритма генерации таблицы, сложность которого зависит, собственно говоря, только линейно от размера этой таблицы. Лично у меня, в общем, не оптимизированный вариант наложения цвета зверски тормозил. SSE снова подняло производительность, я бы сказал до небес. Давайте теперь, когда мы уже примерно представляем себе требуемый алгоритм, взглянем на ещё парочку ассемблерных файликов, занимающихся, собственно, быстрым наложением цветов, начнём, собственно, с renderer_fast.asm:


Код: Delphi
push eax
push ebx
push ecx
push edx

push edi
push esi

movss xmm6, const1

//раскопируваем ее по всему регистру
shufps xmm6, xmm6, 0
movss xmm0, const256
movss xmm6, xmm0
shufpd xmm6, xmm6, 1

mov edx, zskip

mov ebx, cl_gradiento_float
mov edi, line
mov esi, map

mov eax, wy
sal eax, 2  //умножаем на 4
add esi, eax

mov ecx, w
//делим на 4
sar ecx, 2

@for_row:

	movaps xmm0, [esi]

	//по быстрому сначала считаем дробные части
	cvttps2dq xmm1, xmm0
	cvtdq2ps xmm1, xmm1

	//теерь в xmm0 храняться коэффициенты интерполяции
	subps xmm0, xmm1

	//грузим единицы, и вычитаем из них коэффициенты
	movaps xmm2, xmm6
	shufps xmm2, xmm2, 0 //забыли раскопировать лов-парт по всему регистру

	subps xmm2, xmm0
	//меняем местами
	shufpd xmm6, xmm6, 1 

	//в xmm1 - индексы цветов
	//в xmm2 - коэффициенты интерполяции
	//в xmm0 - коэффициенты интерполяции второго порядка
	//в xmm6 - константа, которую надо беречь ака зеницу ока
	//в xmm7 - будем таки накапливать сумму цветов 
	//нам теперь надо просто перемешать нужные цвета

		//PASS-1
		{$I inc/renderer_fast_one_pass.asm}

	shufps xmm0, xmm0, 57
	shufps xmm1, xmm1, 57
	shufps xmm2, xmm2, 57


		//PASS-2
		{$I inc/renderer_fast_one_pass.asm}

	shufps xmm0, xmm0, 57
	shufps xmm1, xmm1, 57
	shufps xmm2, xmm2, 57

		//PASS-3
		{$I inc/renderer_fast_one_pass.asm}

	shufps xmm0, xmm0, 57
	shufps xmm1, xmm1, 57
	shufps xmm2, xmm2, 57
	
		//PASS-4
		{$I inc/renderer_fast_one_pass.asm}

	//в самом конце, происходит сброс в буфер вывода, или линию картинки
	movntps [edi], xmm7 
	//меняем местами (возвращаем в исходное положение)
	shufpd xmm6, xmm6, 1 

	add esi, 16
	add edi, 16

	dec ecx
jnz @for_row

pop esi
pop edi

pop edx
pop ecx
pop ebx
pop eax

Выглядит загадочно и непонятно… Впрочем процент матерных комментариев в коде минимален я старался писать максимально разборчиво, так-что по комментариям можно получить довольно полное представление о том, как именно это работает. Регистр XMM6 тут содержит в себе полезные константы, а именно 1.0 и 256.0, которые, по мере нужды в них, ставятся на первое место командой SHUFPD по очереди. С помощью команды SHUFPS и маски 57 (00`11`10`01), реализуется такой-же циклический сдвиг, как и в предыдущем примере. Для округления в нижнюю сторону используются в связке CVTTPS2DQ и CVTDQ2PS, преобразующие 4 штуки float в 4 long'а, и наоборот, соответственно, с отсечением дробной части, то-бишь округлением вниз. Указатель cl_gradiento_float показывает на выровненную в памяти, удобоваримую таблицу наших цветов, содержащих, собственно, уже готовые четвёрки флоатов, которые нам остаётся только загружать.


Непосредственно процесс смешивания цветов вынесен в отдельный файл, и подключен, как выше видно, инлайн 4 раза — такой подход давал небольшой выигрыш перед циклом, за счет экономии на командах и, возможно, суперскалярности. Давайте мы, кстати на него взглянем, он несколько менее таинственен (это прямое включение из файла renderer_fast_one_pass.asm):


Код: Delphi
cvtss2si eax, xmm1

sub eax, edx //вычитаем zskip
//проверяем на ноль
cmp eax, 0


//jns +02
db $79, $02


	xor eax, eax

//трудно догадаться, но это была такая проверка на eax меньше нуля

sal eax, 4 //умножаем на 16

//берем цветовые констаны 
movaps xmm3, [ebx + eax]
add eax, 16
movaps xmm4, [ebx + eax]

//умножаем их на соответствующие коефиценты
movss xmm5, xmm0
shufps xmm5, xmm5, 0

mulps xmm4, xmm5

movss xmm5, xmm2
shufps xmm5, xmm5, 0

mulps xmm3, xmm5

//складываем (это у нас была линейная интерполяция)
addps xmm3, xmm4

//теперь весь xmm3 надо умножить на 255
movss xmm4, xmm6
shufps xmm4, xmm4, 0 //грузим первое число

//умножаем на 256
mulps xmm3, xmm4

//дальше надо все это запаковать в 4 8-мибиттных числа
cvttps2dq xmm3, xmm3
movaps xmm4, xmm3

psrldq xmm4, 3
por xmm3, xmm4

psrldq xmm4, 3
por xmm3, xmm4

psrldq xmm4, 3
por xmm3, xmm4
//теперь в младшей части xmm3 лежит код цвета, который надо сбросить в регистр-акумулятор

movss xmm7, xmm3

shufps xmm7, xmm7, 57 //00,11,10,01 (сдвиг в обратную сторону)

Тут была съедена интересной породы собака: дело в том, что когда я попытался в файле поставить метку, я слегка не учёл, что подключается он 4 раза, и поэтому Delphi справедливо отметила, что метка эта дублируется. Пришлось обходиться без метки… Впрочем, писать на ассемблере в опкодах — занятие местами захватывающее и интересное :). А вообще, если сказать коротко, то здесь берутся 2 соседних элемента из таблицы цветов, над ними по вычисленным ранее коэффициентам производиться линейная интерполяция, и результат, с помощью MMX, битовых сдвигов, округления, и ядрёной матери, преобразуется в четверку байтов, и пихается в регистр-накопитель, который за этим циклически сдвигается. Всё просто как в швейцарских часах :). Получается, в итоге, что за раз, весь этот алгоритм, какбэ обрабатывает сразу 4 пикселя, ну, за 4, собственно, прохода.


По конкретнее тут можно рассмотреть интересные MMX-подобные команды PSRLDQ, и POR. Первая занимается побайтовым сдвигом регистра на указанное количество байтов, а вторая — ни что иное, как побитовый логический OR. Если приглядеться, то видно, что с помощью этих команд производиться упаковка единственных значащих байтов DWORD чисел, которые остаются в регистрах после того как срабатывает CVTTPS2DQ, в одно-единственное DWORD число. Более умного способа это сделать, я, честно говоря, не нашёл. Или в SSE нету команды, которая переделывает 4 флоата в 4 байта, или я плохо искал. Впрочем, работает оно и так, и работает, на самом деле, достаточно быстро. Кстати использование здесь именно CVTTPS2DQ не совсем оправдано. По хорошему, нужно правильно настроить округление в регистрах состояния, и использовать CVTPS2DQ (с одной T), но мне, видимо, было лень. Ну ничего, покажем, как это сделать, в следующий раз. :)


Множество Мандельброта. Зимний стиль
Рис. №889. Множество Мандельброта. Зимний стиль

Развязка: срываем покровы.


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


Серые кружева
Рис. №890. Серые кружева

Тут уже всё зависит от фантазии и вдохновения…


Красно-чёрный стиль
Рис. №891. Красно-чёрный стиль

Все эти образцы выполнены уже с косинусной интерполяцией и 16-кратным сглаживанием.


Ещё один эксперимент с цветом
Рис. №892. Ещё один эксперимент с цветом

Расскажем теперь напоследок о самой программе, которую, кстати, наконец-то можно скачать вместе с исходниками. Во первых, программа может не заработать, если у вас отсутствует поддержка SSE4 процессором. Но, думаю, наши времена уже давно миновали каменный век, и SSE4, являющееся в общем-то уже довольно старым, должно присутствовать практически повсеместно. Ну, мы на это по крайней мере надеемся. Далее, небольшой экскурс, как программой пользоваться. У неё есть 2 режима: произвольная навигация по фракталу, и работа с буфером-снимком, который можно сделать произвольного размера, глубины, и кратности сглаживания. Навигация производиться мышью и колёсиком, рекомендую ставить уровень прорисовки 4 … 8 на главной вкладке, иначе прорисовка будет очень медленной. Сглаживание выше 16 ставить не рекомендую — скорее всего, не хватит памяти.


Вкладка «Color»
Рис. №893. Вкладка «Color»

Вкладка Color позволяет настраивать таблицу цветов с помощью опорных точек, которые там мы видим в списке. Оные можно удалять/добавлять 4-мя кнопочками под списком, выбрав конкретный узел — менять его цвет, и длину, двигая ползунки ниже. Пусть вас не смущает, что ползунков цвета там 4 — последний был зарезервирован для альфа-канала, и в целом не используется. Таблицу узлов можно копировать в буфер обмена и обратно с помощью кнопок IN и OUT. Чуть ниже можно выбрать тип интерполяции, и установить галочку Cycle, которая позволяет зацикливать вашу палитру по фракталу, указанное в поле ниже число раз. Скриншот можно сделать с помощью кнопки Snap. Ещё там есть важный ползунок Skip, рекомендуется нажимать кнопочку Auto рядом с ним, если вдруг ваш фрактал выглядит как-то странно. Настраивает этот ползунок минимальную высоту на карте, которую нужно вычитать из всех высот при наложении цвета, программа умеет вычислять её автоматически.


Вкладка «Render»
Рис. №894. Вкладка «Render»

На вкладке Render представлено всё, что может понадобиться для работы с буферами-снимками. Алгоритм таков: на вкладке Main путешествуем по фракталу, выбираем понравившуюся область, открываем Render, и делаем снимок выбранной области, в желаемом разрешении и сглаживании, с помощью любой из двух кнопок вверху. Отличие только в том, что правая кнопка устанавливает размеры буфера равные размеру окна программы, а левая берёт размер из выбранной кнопочки в Resolution. Буфер можно сохранить, и потом загрузить, и даже переименовать соответствующими кнопочками. Если во время рендера буфера нажать и держать Esc, то рендер остановиться — полезно, если вам внезапно надоест ждать. Ах да, самое главное! Ползунок Depth задаёт, собственно, глубину прорисовки, его советую ставить на величины 3000 … 6000. Действует этот ползунок только на рендер, в предпросмотре на главной вкладке глубина считается автоматически, на основе длинны цветовой таблицы.


Главная вкладка
Рис. №895. Главная вкладка

Ну и на главной вкладке, как уже писалось, можно по множеству путешествовать. Approximation level задаёт точность прорисовки, рекомендую ставить в пределах 4 – 8, выбирать нужно в зависимости от того, как медленно будет всё рисоваться. Вверху в полях находятся ваши координаты на комплексной плоскости, их можно скопировать соответствующей кнопкой, а можно вставить свои координаты и нажать GO! — программа услужливо перейдёт на указанные координаты. Кнопкой Snap можно тоже сделать скриншот.


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


Код: Delphi
TTraceBufferHeader = packed record
  magic:DWORD;  //в лучших традициях
  width, height, aa_level, tdepth: LongInt;
  re, im, zoom:Double; //лучше это сохранять)
  map:array[0..0]of single;
end;

Сразу за структурой идут непосредственно данные, то-есть наша карта, размером в байтах W ⋅ H ⋅ 4. Никаких сложностей, никакого сжатия, всё проще чем даже формат BMP. Кстати для значения поля magic, отвечающего за идентификацию нашего “формата” используется численная константа 777`777`777, или последовательность байтов 0x71, 0xF2, 0x5B, 0x2E. Впрочем это уже так, для полноты картины.


Так-же пару слов нужно сказать про большую чёрную полосу с кривой на вкладке Render. Я думаю многие догадались, что это ни что иное, как гистограмма карты фрактала, ну, или нечто на нее похожее. Строиться эта гистограмма с помощью так называемого “массива частот”, опять-же с поправкой на то, что пиксели у нас имеют тип float. Высота на графике вычисляется по значениям этого массива по следующей формуле:


\[ H_i = \log^2{\sqrt{x_i}} = {{\log^2{x_i}} \over {4}} \]
Визуализация гистограммы

Здесь Hi означает высоту точки кривой на графике, а xi — количество подсчитанных пикселей с i-тым значением. Формула подобрана экспериментально, на мой взгляд так кривая гистограммы смотрелась наиболее доходчиво. Дело в том, что без использования логарифма, пиксели, которых на карте было большинство, выдавали гигантские пики на кривой, по сравнению с которыми любые мелкие детали становились ничтожно маленькими, и это приводило к тому, что кривая гистограммы вырождалось во что-то уродливое, и совершенно не информативное. Улучшить положение можно было, в принципе, с помощью корней больших степеней, но логарифмы, как мне показалось, выдавали более разборчивый результат. Так их, в общем, я и оставил…


Красно-чёрный
Рис. №429. Ещё один примерчик

Вот мы и завершили рассмотрение нашей великой программы, и пришло время для главного — сейчас мы раскроем её полные исходники, и уже скомпилированный .exe, который должен работать на большинстве Windows. Кстати, почему я её назвал именно Fractalize++, сам не знаю, это было скорее всего первое что пришло в голову :). Хотя два плюса в конце весьма оправданы, это довольно расширенная версия всего того, что я писал для генерации фракталов до этого. И не ругайтесь за то что исходники на Делфи — ну люблю я этот хорошенький, уютненький язычок, и все тут :). И ещё чтобы проект скомпилировался у вас, нужно будет скорее всего правильно настроить пути.


  • Скомпилированная, стабильная версия программы (требует поддержки SSE4+)
  • Исходники программы (Delphi 7);


В следующей статье по этой теме мы рассмотрим некоторые вещи, о которых говорили выше, главная из которых, которая меня на данный момент очень интересует, это, несомненно, распараллеливание процесса вычислений карты, посредством исключительно Windows API (в частности функции WaitForMultipleObjectsEx), и прямом использовании CPU Affinity Mask для гарантированного исполнения потоков на разных процессорах. Может придумаем и что-нибудь ещё :) Ещё несколько образцов, напоследок:


На тему космоса
Рис. №431. На тему космоса

Смелая расцветка
Рис. №436. Смелая расцветка

Bloody-style
Рис. №433. Bloody-style

На этом — до скорых встреч, друзья!

Ключевые слова Программизм, математика, фрактал, Множество Мандельброта, SSE, Оптимизация, ADDSUBPD, Delphi 7, DPPD, интерполяция, гистограмма
QR-код
Интересное
Неприятный сюрприз от Hitach Ultrastar Не долго раздумывая, включил Викторию версии 4.3, и посмотрел SMART. Признаться, не смотря на то что Вика показала зеленый GOOD, смарт меня весьма обескуражил

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