Генерация фракталов с помощью SSE/SSE2 - Алгоритм

ОГЛАВЛЕНИЕ

Алгоритм

Ниже приведен алгоритм для генерации набора Мандельброта в псевдокоде:

комплексные P, Z;
для каждой точки P на комплексной плоскости {
   Z = 0;
   for i = 0 to ITER {
      if abs(Z) > 4.0
         установить цвет точки на palette[i];
         прервать внутренний цикл и перейти к следующей точке;
      Z = Z * Z + P;
   }
   установить цвет точки на черный;
}

P и Z – комплексные переменные. Можно прочитать общие объяснения о комплексных числах в вышеназванных статьях.

Для каждого пикселя P изображения цикл повторяется ITER раз. Выполнение цикла занимает O(X*Y*ITER) времени, где ITER – количество итераций, а размеры изображения равны X x Y пикселей.

abs(Z) – функция, вычисляющая модуль комплексного числа. Представьте круг с радиусом 2.0 и центром в начале координат. Если точка Z находится вне этого круга, abs(Z) будет больше 4.0. В этом случае внутренний цикл должен прерваться, потому что известно, что каждая точка набора Мандельброта лежит внутри круга. И если Z вышло за пределы круга, то точка P не принадлежит набору Мандельброта. Отметим эту точку цветом, показывающим, как быстро она покидает круг при повторении вычисления Z = Z * Z + P. Количество полных итераций i показывает это. Берется i-й цвет из некоторой палитры и продолжается со следующей точки. Палитра – массив цветов. Например, точки, выходящие за пределы круга после первой итерации, отмечаются темно-зеленым цветом; точки, покидающие круг после второй итерации, отмечаются светло-зеленым, и т.д.

Иначе говоря, если цикл выполнился ITER раз, и точка Z вообще не вышла за пределы круга, то точка принадлежит набору Мандельброта. Она окрашивается черным цветом.
Внутренний цикл многократно вычисляет Z<SUB>N+1</SUB> = Z<SUB>N</SUB>2 + P, итерационную формулу для набора Мандельброта. Z инициализируется нулем, но легко понять, что можно установить Z = P, потому что Z<SUB>1</SUB> = Z<SUB>0</SUB>2 + P = 02 + P = P.
Использование палитры – простейший и быстрейший способ окрашивания в наборе Мандельброта, но полученное изображение будет иметь только ITER+1 цветов. Есть методы отрисовки 24-битных изображений. Например, Chaos Pro использует следующую формулу:

color_index = (i - log(log(abs(Z))) * 1.442695 + 2) / 252.0;

Коэффициенты 252.0 и 2.0 уменьшают значение color_index, чтобы оно вмещалось в диапазон 0.1. Затем color_index переводится в 24-битный цвет. Увы, вычисление логарифма – медленная операция, увеличивающая общее время выполнения на 30-50%. Поэтому от 24-битного цвета отказались ради скорости и выбрали ограниченную палитру цветов.

Чтобы перевести алгоритм в код C, надо помнить формулы для умножения комплексных чисел и нахождения их модулей.

void PaintMandelbrotFPU(int w, int h, long* b, double left, double top,
                        double right, double bottom) {
   int x, y, i; double x2, y2, px, py, zx, zy, dx, dy;
   dx = (right - left) / w; // Найти горизонтальное и вертикальное увеличение для каждого пикселя
   dy = (bottom - top) / h;
   py = top;
   for(y = 0; y < h; ++y) {
      px = left;
      for(x = 0; x < w; ++x) {
         zx = px, zy = py;
         for(i = 0; i < ITER; ++i) {
            x2 = zx * zx,   y2 = zy * zy;
            if(x2 + y2 > 4.0)
               break;
            zy = zx * zy * 2 + py; // Вычислить z = z * z + p
            zx = x2 - y2 + px;
         }
         *b++ = a[i];
         px += dx; // Перейти к следующему пикселю
      }
      py += dy; // Перейти к следующей строке пикселей
   }
}

Комплексная переменная p превратилась в px (действительная часть) и py (мнимая часть). То же самое происходит для zx и zy. Квадраты zx и zy используются дважды: при вычислении модуля и при нахождении z * z + p. Поэтому для них было сделано две переменных x2 и y2. Переменные zx и zy инициализируются значениями px и py вместо нулей, потому что пропуск умножения 0*0 экономит время.

Обратите внимание на целые переменные x и y, считающие координаты текущего пикселя, и на пару (px, py), являющуюся координатами на комплексной плоскости. px и py увеличиваются независимо от x и y. Это кардинально повышает производительность, потому что не надо переводить экранные координаты (x, y) в координаты комплексной плоскости (px, py) при каждом повторении цикла. Плата за такое увеличение производительности – потеря точности. Числа с плавающей точкой имеют ограниченную точность, и ошибки накапливаются в ходе многократных сложений: px += dx, py += dy. Это вызывает неприятный эффект: при увеличении масштаба изображения оно будет дергаться. Альтернативным решением могло быть вычисление px = (right - left) * x / w, py = (bottom - top) * y / h на каждой итерации, но оно непрактично из-за большой задержки на деление.

Наконец, a – палитра цветов, а b – битовая карта, генерируемая функцией. Битовая карта b позже переносится на экран с помощью функций DirectDraw или GDI(интерфейс графического устройства). Детали этого процесса неважны. Просто помните, что размер b равен w * h пикселей, и что надо заполнить его цветами из палитры.

Когда внутренний цикл прерывается с помощью break, значение i меньше ITER. А когда он не прерывается, i == ITER после цикла. Благодаря данному свойству был создан массив a[ITER+1], и последнему элементу этого массива был присвоен черный цвет. После цикла a[i] ссылается на черный a[ITER]. Более строгий способ сделать это был бы:

// ...
for(i = 0; i < ITER; ++i) {
    x2 = zx * zx,   y2 = zy * zy;
    if(x2 + y2 > 4.0) {
       *b++ = a[i];
       px += dx;
       goto NEXTPIXEL;
    }
    zy = zx * zy * 2 + py; // Вычислить z = z * z + p
    zx = x2 - y2 + px;
}
*b++ = RGB(0,0,0);
px += dx; // Перейти к следующему пикселю
NEXTPIXEL:
// ...

Алгоритм для генерации набора Джулии схож с алгоритмом для набора Мандельброта, поэтому не показан тут. Главное отличие – введение комплексного параметра (PX, PY) в Джулии.