Зона кода

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

Построение изображений с помощью итерационных функций

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

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

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

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

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

Итак, переходим к теоретической части нашего повествования.

Итерационные функции и случайные последовательности

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

Зададим 2 последовательности, xnn=1 и ynn=1, с помощью следующих рекуррентных формул:

xn=fxn-1,yn-1, n, yn=gxn-1,yn-1, n.

Поясним, что x0 и y0 — это некоторые заранее заданные числа, а f(xy) и g(xy) — это некоторые функции двух переменных, называемые итерационными. Сам процесс вычисления очередного члена той или иной последовательности через такие функции будем называть итерациями, а приведённый выше набор рекуррентных формул — итерационной схемой.

Рекурсивный способ задания последовательностей, скорее всего, хорошо знаком читателю, если он изучал математику в вузе. Несколько необычным может показаться "перекрёстный" способ вычисления членов последовательностей, при котором для вычисления n-го члена каждой из двух последовательностей нужен не только n − 1-й член той же последовательности, но и n − 1-й член другой.

А теперь рассмотрим схему построения членов двух последовательностей, использующую не одну пару итерационных функций, а m пар. Каждая из этих функций будет линейной по обеим переменным, а также будет содержать аддитивную константу. Более конкретно, функции будут иметь вид:

fkx,y=akx+bky+ckgkx,y=dkx+eky+hk, k=0,1,,m-1.

Для каждого n, начиная с 1, будет случайным образом выбираться число от 0 до m − 1, и при вычислении xn и yn в рекуррентных формулах будет использоваться пара итерационных функций, индексы которых равны данному случайному числу. Отметим, что случайные числа, "появляющиеся" перед каждой итерацией, не обязаны быть равновероятными. Однако для разных шагов вероятность появления конкретного фиксированного числа одна и та же.

Давайте теперь сформулируем сказанное на строгом математическом языке. Рассмотрим последовательность дискретных независимых в совокупности случайных величин Tn=1 , распределённых по одному и тому же закону. А именно: каждая случайная величина принимает значения 0, 1, …, m − 1 с соответствующими вероятностями p0p1, …, pm-1.

Теперь последовательности, xnn=1 и ynn=1 зададим с помощью следующей итерационной схемы:

xn=fTnxn-1,yn-1, n, yn=gTnxn-1,yn-1, n.

Как и ранее, x0 и y0 — это некоторые заранее заданные числа.

Таким образом, каждая из последовательностей является случайной, т. е. её члены — это случайные величины. Однако, каждую из этих последовательностей можно "реализовать", т. е. вычислить все её члены (разумеется, таких реализаций будет бесконечно много).

Зададимся главным вопросом данного раздела. А какое же отношение изображения, которые мы собираемся строить, имеют к этой паре случайных последовательностей? Очень простое. Построим реализацию этих двух последовательностей. Для каждого натурального n пару (xnyn) можно рассматривать как координаты точки, заданной в декартовой прямоугольной системе координат на плоскости. Так вот, изображение, соответствующее некоторой паре реализованных последовательностей, представляет собой геометрическое место всех таких точек на плоскости.

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

Добавим, что при построении изображений на компьютере, мы, разумеется, будем выполнять лишь конечное (но достаточно большое) число итераций.

О генерации псевдослучайных чисел

При написании программы мы столкнёмся с необходимостью генерировать псевдослучайные числа, распределённые, вообще говоря, не равномерно, а по заранее заданному закону. В то же самое время, мы будем располагать лишь программным генератором псевдослучайных чисел, равномерно распределённых на промежутке [0, 1]. Как из второго распределения получить первое?

Переведём задачу в математическую плоскость. Пусть имеется непрерывная случайная величина U, распределённая равномерно на отрезке [0, 1]. Зададимся целью построить дискретную случайную величину T как функцию от U, таким образом, чтобы T принимала значения 0, 1, …, m − 1 с соответствующими вероятностями p0p1, …, pm-1.

Решить поставленную задачу весьма просто. Введём в рассмотрение суммы вероятностей

sk=i=0k-1pi, k=0,1,,m-1.

Если верхний предел суммирования по i меньше нижнего, то такую сумму по определению будем полагать равной 0.

Т выразим через U следующим образом:

T=0, если Us0,s1,1, если Us1,s2,2, если Us2,s3,     ,     ,m-1, если Usm-1,1.

Очевидно, случайная величина T распределена по требуемому нами закону. Заметим, что, по сути, Т — это номер промежутка, в который попадает случайная величина U (при условии, что промежутки мы нумеруем числами от 0 до m − 1 в порядке возрастания их левых границ).

С практической точки зрения полученный результат позволяет на каждом шаге итерации в качестве номера итерационных функций брать номер промежутка, в который попадает число, сгенерированное датчиком псевдослучайных чисел, равномерно распределённых на отрезке [0, 1].

А теперь можно переходить к написанию программы.

Структура программы

Программа состоит из файла main.c и файлов, образующих графическую библиотеку pgraph. Содержимое файла main.c начинается со следующей директивы, подключающей графическую библиотеку:

#include "pgraph.h"

Далее в файле содержатся описания глобальных константных переменных и константных массивов. За ними — определения функций get_random_value() и main(). Первая из них генерирует псевдослучайные числа, а вторая выполняет основную работу по построению изображений.

Глобальные константные переменные и константные массивы

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

Ниже приводятся описания данных констант и массивов.

  • n — количество итераций;
  • w — ширина изображения в пикселях;
  • h — высота изображения в пикселях;
  • xc — абсцисса начала новой системы координат в старой системе;
  • yc — ордината начала новой системы координат в старой системе;
  • l — длина в пикселях отрезка, параллельного одной из осей координат, имеющего в новой системе координат единичную длину;
  • m — количество пар итерационных функций, т. е. число m;
  • s — одномерный массив размера m, содержащий суммы вероятностей случайных величин Tn (k-й элемент массива содержит sk);
  • f — двухмерный массив, состоящий из m "строк" и 3-х "столбцов", содержащий константы, задействованные в функциях fk(xy) (элементы массива с индексами (k, 0), (k, 1), (k, 2) содержат числа ak, bk, ck соответственно, где 0 ≤ k ≤ m − 1);
  • g — двухмерный массив, состоящий из m "строк" и 3-х "столбцов", содержащий константы, задействованные в функциях gk(xy) (элементы массива с индексами (k, 0), (k, 1), (k, 2) содержат числа dk, ek, hk соответственно, где 0 ≤ k ≤ m − 1).

Все переменные имеют тип int, а базовым типом всех массивов является double.

Поясним, что под "старой" системой координат подразумевается та, которая определена в библиотеке pgraph. Построения всех изображений будут вестись в новой системе, полученной из старой параллельным переносом (сдвиги по осям абсцисс и ординат равны соответственно xc и yc) и "сжатием" в l раз. Таким образом, точка, имеющая в новой системе координаты (xy), в старой будет иметь координаты (x l + xcy l + yc). Излишне, думаю, пояснять, что за хранение чисел xc, yc и l ответственны константные переменные xc, yc и l соответственно.

Для хранения чисел x0 и y0 переменные не выделяются, поскольку во всех случаях построения изображений в качестве этих чисел берутся нули.

Генерация псевдослучайных чисел: функция get_random_value()

Функция get_random_value() при каждом обращении к ней генерирует псевдослучайное целое число в диапазоне от 0 до m − 1 в соответствии с описанной ранее схемой. Вот код этой функции:

1.int get_random_value()
2.{
3.    double r = (double) rand() / RAND_MAX;
4.    int c = 1;
5.    while (s[c] < r && ++c < m)
6.        ;
7.    return c - 1;
8.}

Получаем с помощью стандартной библиотечной функции rand() псевдослучайное число в диапазоне от 0 до значения макроса RAND_MAX, делим полученный результат на это значение и присваиваем частное переменной r (стр. 3). Теперь в r хранится число, принадлежащее отрезку [0, 1]. Его приближённо можно считать значением случайной величины, равномерно распределённой на этом отрезке.

Поясним, что значение макроса RAND_MAX, в нашем случае (т. е. в случае использования компилятора MinGW64 версии 4.9.2 для 64-битных систем) равно 32767.

Теперь, с помощью линейного поиска, задействующего цикл while, ищем индекс наибольшего элемента массива s, не превосходящего значение r, увеличенный на единицу, и сохраняем его в переменной c (см. стр. 4-6). Отметим, что в случае, если значение r — нулевое, цикл не выполняется ни разу, а переменная с сохраняет единичное значение (см. стр. 4).

Теперь значение c, уменьшенное на единицу, которое мы возвращаем в строке 7, является номером того промежутка из набора [s0s1], (s1s2], (s2s3], …, (sm-1, 1], в который попадает значение r (за подробностями обращайтесь к разделу "О генерации псевдослучайных чисел").

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

Генерация изображения: функция main()

А вот и код функции main():

 1.int main()
 2.{
 3.    image *img = create_image(w, h);
 4.    double x = 0, y = 0;
 5.    for (int i = 0; i < n; i++)
 6.    {
 7.        int r = get_random_value();
 8.        double x1 = f[r][0] * x + f[r][1] * y + f[r][2];
 9.        double y1 = g[r][0] * x + g[r][1] * y + g[r][2];
10.        x = x1;
11.        y = y1;
12.        set_color(img, round(x * l) + xc, round(y * l) + yc, BLACK);
13.    }
14.    save_to_file(img, "out.bmp");
15.    free(img);
16.    return 0;
17.}

Создаём изображение с заданными размерами (стр. 3). Выделяем память под переменные x и y, в которых будут храниться текущие члены последовательностей, и инициализируем их нулями (стр. 4). Напомню, что в качестве чисел x0 и y0, участвующих в вычислении первых членов каждой из последовательностей, берутся нули.

Вычисляем в цикле for первые n членов каждой последовательности (стр. 5-13). Получаем сначала псевдослучайное число и записываем его в r (стр. 7). Далее вычисляем текущие значения членов обеих последовательностей, помещая их во временные переменные x1 и y1 (стр. 8, 9). При вычислении используем константы, фигурирующие в итерационных функциях и хранящиеся в массивах f и g. Выбор той или иной пары наборов коэффициентов (а значит, пары итерационных функций) зависит от значения r, использующегося в качестве первых индексов участвующих в вычислениях элементов массивов.

Переписываем вычисленные текущие значения в переменные x и y (стр. 10, 11). Координаты точки, содержащиеся в этих переменных, переводим в координаты исходной системы координат, округляем до целых и наносим точку с результирующими координатами на изображение чёрным цветом (стр. 12).

По завершении цикла сохраняем сформированное изображение в файле "out.bmp" (стр. 14) и освобождаем занимаемую изображением память (стр. 15). На этом работа функции завершается.

Построение изображения треугольника Серпиньского

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

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

В книге Седжвика и других авторов предлагается следующий способ построения изображения треугольника Серпиньского. Рассмотрим 3 точки на плоскости, являющиеся вершинами равностороннего треугольника, например, точки с координатами 0,0, 0,1, 1/2,3/2 в декартовой прямоугольной системе координат. Выбираем наугад (с равными вероятностями) одну из трёх вершин треугольника и строим точку, делящую отрезок, соединяющий вершину с координатами 0,0 и выбранную наугад вершину, пополам. Это первая точка нашего изображения.

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

Повторяем данный процесс многократно (авторы книги предлагают строить в общей сложности 100000 точек), в результате чего получаем изображение треугольника Серпиньского.

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

Нам потребуются 3 пары итерационных функций. Их индексы 0, 1, 2 должны выбираться с вероятностями 1/3, 1/3, 1/3 соответственно. Сами итерационные функции приведены ниже.

f0x,y=1/2x,g0x,y=1/2y,f1x,y=1/2x+1/2,g1x,y=1/2y,f2x,y=1/2x+1/4,g2x,y=1/2y+3/4.

Теперь давайте вставим в нашу программу описания глобальных константных переменных и константных массивов, соответствующие данным вероятностям и данным итерационным функциям. Но для начала определим макрос TRIANGLE, поместив в файл main.с после инструкции #include следующую инструкцию

#define TRIANGLE

После инструкции вставляем в файл следующий код:

//Треугольник Серпиньского
#ifdef TRIANGLE
const int n = 100000;                                  //количество итераций
const int w = 620, h = 550;                            //размеры изображения
const int xc = 10, yc = 10;                            //координаты начала новой системы координат в старой
const int l = 600;                                     //коэффициент сжатия
const int m = 3;                                       //количество пар итерационных функций
const double s[] = {0, 0.3333333, 0.6666667};          //массив сумм вероятностей
const double f[][3] = {{0.5, 0.0, 0.0},                //массив коэффициентов для функций f(x,y),
                       {0.5, 0.0, 0.5},                //задействованных для вычислений x
                       {0.5, 0.0, 0.25}};
const double g[][3] = {{0.0, 0.5, 0.0},                //массив коэффициентов для функций g(x,y),
                       {0.0, 0.5, 0.0},                //задействованных для вычислений y
                       {0.0, 0.5, 0.4330127}};
#endif

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

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

Изображение треугольника Серпиньского

Изображение треугольника Серпиньского

Построение изображения папоротника Барнсли

Следующее изображение, построение которого описывается в книге Седжвика и других, — это изображение папоротника Барнсли. Теперь нам уже потребуются 4 пары итерационных функций. Их индексы 0, 1, 2, 3 будут выбираться с вероятностями 0,01, 0,85, 0,07, 0,07 соответственно. А вот и сами итерационные функции:

f0x,y=0,5,g0x,y=0,16y,f1x,y=0,85x+0,04y+0,075,g1x,y=-0,04x+0,85y+0,18,f2x,y=0,2x-0,26y+0,4,g2x,y=0,23x+0,22y+0,045,f3x,y=-0,15x+0,28y+0,575,g3x,y=0,26x+0,24y-0,086.

Вносим теперь изменения в программу. Инструкцию #define заменяем инструкцией

#define FERN

А после #ifdef-блока помещаем следующий фрагмент кода:

//Папоротник Барнсли
#ifdef FERN
const int n = 100000;
const int w = 620, h = 620;
const int xc = 0, yc = 10;
const int l = 600;
const int m = 4;
const double s[] = {0, 0.01, 0.86, 0.93};
const double f[][3] = {{0.0, 0.0, 0.5},
                       {0.85, 0.04, 0.075},
                       {0.2, -0.26, 0.4},
                       {-0.15, 0.28, 0.575}};
const double g[][3] = {{0.0, 0.16, 0.0},
                       {-0.04, 0.85, 0.18},
                       {0.23, 0.22, 0.045},
                       {0.26, 0.24, -0.086}};
#endif

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

Изображение папоротника Барнсли

Изображение папоротника Барнсли

Построение изображения дерева

Теперь построим то, что в книге Седжвика и других авторов называется "деревом", хотя то, что оказывается изображённым, скорее, похоже на набор деревьев различных размеров. На этот раз в итерационном процессе будут участвовать 6 пар итерационных функций. Их индексы 0, 1, 2, 3, 4, 5 будут выбираться с вероятностями 0,1, 0,1, 0,2, 0,2, 0,2, 0,2 соответственно. Вот эти функции:

f0x,y=0,55,g0x,y=0,6y,f1x,y=-0,05x+0,525,g1x,y=-0,5x+0,75,f2x,y=0,46x-0,15y+0,27,g2x,y=0,39x+0,38y+0,105,f3x,y=0,47x-0,15y+0,265,g3x,y=0,17x+0,42y+0,465,f4x,y=0,43x+0,26y+0,29,g4x,y=-0,25x+0,45y+0,625,f5x,y=0,42x+0,26y+0,29,g5x,y=-0,35x+0,31y+0,525.

Заменяем инструкцию #define инструкцией

#define TREE

За последним #ifdef-блоком вставляем следующий код:

//Дерево
#ifdef TREE
const int n = 100000;
const int w = 620, h = 620;
const int xc = 0, yc = 10;
const int l = 600;
const int m = 6;
const double s[] = {0, 0.1, 0.2, 0.4, 0.6, 0.8};
const double f[][3] = {{0.0, 0.0, 0.55},
                       {-0.05, 0.0, 0.525},
                       {0.46, -0.15, 0.27},
                       {0.47, -0.15, 0.265},
                       {0.43, 0.26, 0.29},
                       {0.42, 0.26, 0.29}};
const double g[][3] = {{0.0, 0.6, 0.0},
                       {-0.5, 0.0, 0.75},
                       {0.39, 0.38, 0.105},
                       {0.17, 0.42, 0.465},
                       {-0.25, 0.45, 0.625},
                       {-0.35, 0.31, 0.525}};
#endif

Результат работы скомпилированной программы — это изображение, приведённое ниже:

Изображение дерева

Изображение дерева

Изображение коралла

Последнее изображение, которое мы построим, руководствуясь книгой Седжвика, — это изображение коралла. Нам потребуются 3 пары итерационных функций. Их индексы 0, 1, 2 будут выбираться с вероятностями 0,4, 0,15, 0,45 соответственно. Итерационные функции приведены ниже.

f0x,y=0,3077x-0,5315y+0,8863,g0x,y=-0,4615x-0,2937y+1,0962,f1x,y=0,3077x-0,0769y+0,2166,g1x,y=0,1538x-0,4476y+0,3384,f2x,y=0,5455y+0,0106,g2x,y=0,6923x-0,1958y+0,3808.

Заменяем инструкцию #define инструкцией

#define CORAL

За последним #ifdef-блоком вставляем новый блок:

//Коралл
#ifdef CORAL
const int n = 100000;
const int w = 620, h = 620;
const int xc = 10, yc = 10;
const int l = 600;
const int m = 3;
const double s[] = {0, 0.4, 0.55};
const double f[][3] = {{0.3077, -0.5315, 0.8863},
                       {0.3077, -0.0769, 0.2166},
                       {0.0, 0.5455, 0.0106}};
const double g[][3] = {{-0.4615, -0.2937, 1.0962},
                       {0.1538, -0.4476, 0.3384},
                       {0.6923, -0.1958, 0.3808}};
#endif

Вот какое изображение получаем в результате компиляции и выполнения программы:

Изображение коралла

Изображение коралла

Заключение

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

Отмечу, что, судя по изображениям, треугольник Серпиньского и папоротник Барнсли являются фракталами. Скорее всего, то же самое можно сказать про "дерево" и "коралл", но их фрактальная природа, пожалуй, чуть менее очевидна.

По приведённой ниже ссылке, как всегда, можно скачать исходный код рассмотренной в статье программы. В файле main.c имеются четыре инструкции #define, каждая из которых соответствует одному из четырёх изображений. Три из них закомментированы. Ясно, что для того, чтобы перейти от одного изображения к другому, требуется закомментировать незакомментированную инструкцию и раскомментировать одну из закомментированных. Ну, Вы поняли...улыбка

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

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