Зона кода

А давайте немного попрограммируем...

Визуализация множества Мандельброта

C99Графические программы

В этой статье мы рассмотрим процесс визуализации множества Мандельброта. Идея статьи пришла мне в голову достаточно давно, но некий толчок к её написанию я получил после того, как полистал книгу Седжвика, Уэйна и Дондеро "Программирование на языке Python".

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

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

Возвращаемся к теме статьи. Алгоритм визуализации множества Мандельброта я взял из книги Седжвика и реализовал на языке C99. Кроме того, я последовал предложению автора "раскрасить" изображение множества Мандельброта. В программах, которые будут рассмотрены в статье, используется графическая библиотека pgraph. Предполагается, что читатель знаком с нею хотя бы в общих чертах.

Для понимания нижеизложенного материала также пригодятся элементарные знания теории функций комплексного переменного и теории последовательностей.

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

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

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

Пред тем, как переходить к визуализации, дадим определение множества Мандельброта.

Рассмотрим, вначале, произвольное комплексное число z0 = x0 + iy0. Этому числу соответствует точка на комплексной плоскости с координатами (x0y0). В дальнейшем будем отождествлять комплексные числа с соответствующими им точками комплексной плоскости. Комплексную последовательность

z0, z1, z2, …, где zk = (zk-1)2 + z0, k ∈ ℕ,

назовём, последовательностью Мандельброта для точки z0. Ясно, что каждая такая последовательность является либо ограниченной (по модулю, разумеется), либо неограниченной. Множеством Мандельброта называется множество всех таких точек комплексной плоскости, для которых последовательности Мандельброта ограничены.

Простейшая визуализация множества Мандельброта заключается в следующем. Точкам некоторой прямоугольной области комплексной плоскости естественным образом ставятся в соответствие пиксели холста для рисования. Пиксели, соответствующие точкам плоскости, принадлежащим множеству Мандельброта, закрашиваются чёрным цветом, а остальные — белым. Кстати, известно, что точки множества Мандельброта более или менее "помещаются" в квадрат с центром в точке (-1/2, 0) со стороной 2.

Для построения такой визуализации, очевидно, требуется алгоритм, устанавливающий принадлежность или непринадлежность той или иной точки комплексной плоскости множеству Мандельброта. Проблема заключается в том, что универсального алгоритма, "работающего" для любой точки, не существует. Однако известно, что последовательность Мандельброта, модуль некоторого члена которой превышает 2, стремится по модулю к бесконечности, т. е. не ограничена. Это означает, что начальный член этой последовательности не принадлежит множеству Мандельброта.

Это утверждение позволяет нам использовать алгоритм определения "условной" принадлежности точки данному множеству. Заключается он в вычислении фиксированного, заранее выбранного числа первых членов последовательности Мандельброта для проверяемой точки. Если модуль ни одного из вычисленных членов не превосходит числа 2, то испытуемая точка считается принадлежащей множеству Мандельброта. В противном случае точка считается этому множеству не принадлежащей.

Авторы упомянутой ранее книги предлагают ограничиться нахождением первых 255 членов последовательности (включая саму проверяемую точку). Если все 255 чисел оказываются по модулю не превосходящими 2, то соответствующие им пиксели окрашиваются в чёрный цвет. Если же начальный член (т. е. сама проверяемая точка) превосходит по модулю 2, то соответствующий ей пиксель закрашивается в белый цвет.

Если же установлена непринадлежность точки z0 множеству Мандельброта по причине того, что для некоторого минимального k такого, что 1 ≤ k ≤ 254, выполняется неравенство |zk| > 2, то в этом случае логично окрашивать соответствующий точке пиксель также в белый цвет. Однако авторы книги предлагают иной подход: закрашивать пиксель оттенком серого цвета. При этом цвет закраски должен быть тем ближе к белому, чем меньше k и тем ближе к чёрному, чем больше k. Более конкретно: интенсивность серого цвета предлагается вычислять как разность 255 − k.

В данном случае имеется в виду интенсивность, которая совпадает со значениями каждой из трёх компонент при сохранении цветов в RGB-формате. Напомним, что каждому оттенку серого соответствует именно совпадение компонент.

Таким образом, наша "картинка" не будет содержать лишь два цвета: черный и белый. Но будет включать в себя и точки серого оттенка, тем более тёмные, чем "позже" члены последовательностей Мандельброта для этих точек превосходят по модулю число 2.

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

Структура программы. Директивы препроцессора

Программа будет состоять из главного файла main.c, а также трёх файлов, образующих библиотеку pgraph. К main.c будет подключён заголовочный файл библиотеки pgraph.h с помощью директивы #include. Вот как выглядят директивы препроцессора, с которых начинается файл main.c:

#include <complex.h>
#include "pgraph.h"

#define W 681
#define H 681
#define X0 (3*W/4)
#define Y0 (H/2)
#define L 320.0

Как видим, кроме praph.h, к нашему файлу подключается также заголовочный файл стандартной библиотеки для работы с комплексными числами complex.h (о поддержке языком C99 комплексных чисел (в рамках использования компилятора MinGW64) можно почитать здесь).

Далее с помощью директивы #define определяются 5 макросов, имеющих следующий смысл:

  • W — ширина изображения;
  • H — высота изображения;
  • X0 — абсцисса пикселя в исходной системе координат, в который помещено начало новой системы координат;
  • Y0 — ордината пикселя в исходной системе координат, в который помещено начало новой системы координат;
  • L — длина отрезка, имеющего в новой системе координат единичную длину, выраженная в пикселях.

Поясним, что под исходной системой координат понимается та, которая используется в библиотеке pgraph и позволяет обращаться к пикселям изображения по координатам. А новая система координат — это та, которая вводится нами на комплексной плоскости для построения в ней изображения множества Мандельброта.

Легко заметить, что пиксель, имеющий координаты (i, j) в исходной системе координат будет соответствовать точке с координатами

((i - X0) / L, (j - Y0) / L)

в новой системе координат.

Как мы видим (см. код), наше изображение будет квадратным, отображающим область комплексной плоскости размером примерно 2,1x2,1, а начало координат равноудалено от "нижней" и "верхней" границ квадрата и примерно в 2 раза ближе к "правой" его границе, чем к "левой". Выбор параметров, приводящий к описанному положению вещей, обусловлен упоминавшейся ранее особенностью "локализации" точек множества Мандельброта.

В файле main.c будут содержаться определения двух функций: get_gray_color() и main(). Рассмотрим каждую из них.

Вычисление интенсивности серого цвета пикселя: функция get_gray_color()

Ниже приведён код функции get_gray_color().

 1.uchar get_gray_color(double complex z0)
 2.{
 3.    double complex z = z0;
 4.    for (uchar gray = 255; gray; gray--)
 5.    {
 6.        if (cabs(z) > 2)
 7.            return gray;
 8.        z = z * z + z0;
 9.    }
10.    return 0;
11.}

Функция принимает в качестве параметра точку комплексной плоскости и возвращает "цвет", которым должен быть раскрашен пиксель, соответствующий данной точке. Слово "цвет" помещено в кавычки, т. к., если быть точным, то пиксель раскрашивается в один из оттенков серого цвета (включая белый и чёрный цвета), и данная функция возвращает, на самом деле, интенсивность серого цвета в диапазоне от 0 до 255. Поясним, что тип uchar определён в файле pgraph.h и эквивалентен типу unsigned char.

Алгоритм нахождения интенсивности уже был описан ранее. Текущее значение интенсивности содержится в переменной gray (стр. 4). В ходе цикла for (см. стр. 4-9) gray пробегает значения интенсивности от 255 до 1; также в ходе цикла вычисляются значения последовательности Мандельброта для исходной точки (см. стр. 8). Для каждого значения интенсивности проверяем, превосходит ли текущая точка последовательности Мандельброта число 2. Если превосходит, то прерываем цикл и возвращаем текущую интенсивность (стр. 6-7).

Если же модуль ни одного из вычисленных первых 255 членов последовательности не достигает и не превышает 2, то считаем, что начальная точка принадлежит множеству и возвращаем нулевую интенсивность (стр. 10), соответствующую чёрному цвету.

Построение изображения множества Мандельброта: функция main()

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

 1.int main()
 2.{
 3.    image *img = create_image(W, H);
 4.    for (int i = 0; i < W; i++)
 5.        for (int j = 0; j < H; j++)
 6.        {
 7.            double x = (i - X0) / L;
 8.            double y = (j - Y0) / L;
 9.            double complex z = x + I * y;
10.            uchar gray = get_gray_color(z);
11.            if (gray != 255)
12.                set_color(img, i, j, (color) {gray, gray, gray});
13.        }
14.    save_to_file(img, "out.bmp");
15.    free(img);
16.    return 0;
17.}

Создаём изображение (стр. 3). Перебираем в двух циклах for (см. стр. 4, 5) все его пиксели. Для каждого пикселя вычисляем координаты точки на комплексной плоскости, которой он соответствует (стр. 7, 8). Формируем из координат точки комплексное число (стр. 9). Передаём это число в функцию get_gray_color() и получаем "цвет" текущего пикселя (стр. 10). Если "цвет" пикселя отличен от белого, то закрашиваем пиксель данным цветом (стр. 11-12) (изначально холст для рисования белый, поэтому нет необходимости в окрашивании его пикселей в белый цвет).

Записываем изображение в файл "out.bmp" (стр. 14) и удаляем изображение (стр. 15).

Результат работы программы

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

Чёрно-белое изображение множества Мандельброта

Чёрно-белое изображение множества Мандельброта

Раскрашиваем изображение множества Мандельброта

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

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

Очевидно, что изменять функцию get_gray_color() не требуется. А вот функцию main() нужно переделать. Вот её новая, "цветная" версия:

 1.int main()
 2.{
 3.    uchar red[] = {0, 36, 73, 109, 146, 182, 219, 255};
 4.    uchar *green = red;
 5.    uchar blue[] = {0, 85, 170, 255};
 6.    color colors[256];
 7.    int m = 0;
 8.    for (int i = 0; i < 8; i++)
 9.        for (int j = 0; j < 8; j++)
10.            for (int k = 0; k < 4; k++)
11.                colors[m++] = (color) {red[i], green[j], blue[k]};
12.    image *img = create_image(W, H);
13.    for (int i = 0; i < W; i++)
14.        for (int j = 0; j < H; j++)
15.        {
16.            double x = (i - X0) / L;
17.            double y = (j - Y0) / L;
18.            double complex z = x + I * y;
19.            uchar num = get_gray_color(z);
20.            if (num != 255)
21.                set_color(img, i, j, colors[num]);
22.        }
23.    save_to_file(img, "out.bmp");
24.    free(img);
25.    return 0;
26.}

Список предопределённых цветов, о котором шла речь выше, формируется в виде массива colors в строках 3-11. Принцип формирования достаточно простой: выбираются по 8 оттенков красного и зелёного цветов, а также 4 оттенка синего. Далее эти оттенки смешиваются, в результате чего получаются 256 цветов, коими и заполняется массив colors посредством трёх вложенных циклов (см. стр. 7-11). Интенсивности каждого из трёх цветов примерно равномерно распределяются по спектру интенсивностей данного цвета (см. стр. 3-5).

Теперь значение, возвращаемое функцией get_gray_color(), интерпретируем как индекс элемента массива colors, содержащего цвет, которым нужно раскрасить текущий пиксель (см. стр. 19-21). Заметим, что значениям 0 и 1, возвращаемым функцией, как и в предыдущем ("чёрно-белом") случае, соответствуют чёрный и белый цвета пикселя соответственно.

В результате выполнения новой версии программы получаем следующее изображение:

Цветное изображение множества Мандельброта

Цветное изображение множества Мандельброта

Ясно, что существует огромное количество способов раскраски множества Мандельброта. Для получения некоторых из них можно, например, изменить порядок следования выражений внутри фигурных скобок, отвечающих за значения полей элементов массива colors (см. стр. 11). Существуют, очевидно, 6 вариантов расположения выражений:

1.                colors[m++] = (color) {red[i], green[j], blue[k]};
2.                colors[m++] = (color) {red[i], blue[k], green[j]};
3.                colors[m++] = (color) {blue[k], red[i], green[j]};
4.                colors[m++] = (color) {green[j], blue[k], red[i]};
5.                colors[m++] = (color) {blue[k], green[j], red[i]};
6.                colors[m++] = (color) {green[j], red[i], blue[k]};

Заменяя 11-ю строку функции main() одной из приведённых выше строк (за исключением 1-й), получаем ещё 5 вариантов раскраски. Изображения множества Мандельброта, раскрашенного всеми 6-ю способами, соответствующими 6-ю приведённым выше строкам, изображены ниже.

6 изображений множества Мандельброта

6 изображений множества Мандельброта (щёлкните для увеличения)

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

#define W 1281
#define H 1281
#define X0 (3*W/4)
#define Y0 (H/2)
#define L 640.0

В качестве 11-й строки (см. код функции main()) будем использовать следующую:

                colors[m++] = (color) {red[i], blue[k], green[j]};

После запуска и компиляции программы получаем изображение, которое можно посмотреть по ссылке (откроется в новом окне; исходный файл преобразован в формат JPG; размер JPG-файла составляет 475 КБ). Это же изображение можно получить, щёлкнув по приведённой ниже картинке, представляющей уменьшенную его версию.

Ещё одно изображение множества Мандельброта

Ещё одно изображение множества Мандельброта (щёлкните для увеличения)

И последнее. Из алгоритма раскраски множества Мандельброта следует, что раскраске подвергаются лишь пиксели, соответствующие точкам комплексной плоскости, удалённым от начала координат не более, чем на 2, т. е. лежащими внутри круга радиуса 2 с центром в начале координат или на его границе. Давайте рассмотрим этот круг полностью. Для этого изменим значения макросов, определённых в файле main.c следующим образом:

#define W 1281
#define H 1281
#define X0 (W/2)
#define Y0 (H/2)
#define L 320.0

После запуска и компиляции программы получаем изображение, уменьшенная версия которого приведена ниже (оригинал изображения размером 237 КБ, переведённый в формат JPG, можно получить, щёлкнув по картинке или по ссылке, которая откроется в новом окне).

Множество Мандельброта с высоты птичьего полёта

Множество Мандельброта с высоты птичьего полёта (щёлкните для увеличения)

Заключение

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

Скачать исходный код