Monthly Archives: Апрель 2011

AV-1: SortVis

скриншот / визуализатор сортировок

Введение

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

В качестве основы визуализатора взят код из BY-0.

Исходный код: codepad. Архив откомпилированный exe + freeglut.dll можно скачать по ссылке (для запуска exe должны быть установлены библиотеки Microsoft Visual C++ 2010 redistributable).

Принцип работы

Рассматривается сортировка массива чисел. На экране «число» представляется как прямоугольник определённой длины, а процесс сортировки выглядит как последовательность сравнений (зелёный цвет) и присваиваний (красный цвет), упорядочивающая столбец прямоугольников по убыванию (если смотреть сверху-вниз).

Итак, рассматриваем только те алгоритмы, которым достаточно сравнивать (оператор < ) и присваивать (оператор = ) объекты сортируемой последовательности. Чтобы «записать» последовательность вызова указанных операторов вместо чисел сортируются объекты специального класса — «зонда» Probe, который содержит как числовое значение, так и указатель на очередь событий, куда перегруженные operator< и operator= складывают соответствующие события. Такой подход позволяет визуализировать подходящий алгоритм (например, std::sort), не внося в него изменений.

Для каждого столбца создаётся свой массив и своя очередь событий, используется свой алгоритм сортировки. После применения всех сортировок очереди событий «раскручиваются» с определённой временной задержкой (нажмите «+», чтобы ускорить, «-«, чтобы замедлить), применяя к прямоугольникам на экране действия, ранее произведённые с «зондами».

Детали реализации

Прямоугольник (class Rect) является визуализацией числа, задаётся двумя углами (левый нижний с координатами (x0, y0) и правый верхний с координатами (x1, y1)) и цветом.

Массив прямоугольников (class RectArray), также «столбец», является визуализацией сортируемого массива. Имеет две функции-члена, отвечающих за визуализацию событий:

void assign ( unsigned index, float value );

соответствует присваиванию элементу массива с индексом index значения value,  прямоугольник с этим индексом изменит длину и будет помечен красным цветом;

void compare ( unsigned i1, unsigned i2 );

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

Событие (struct IEvent) — абстрактный базовый класс, наследники которого должны определять метод void apply(RectArray&), вызов которого с передачей соответствующего массива прямоугольников и будет приводить к применению события и его последующей визуализации.

Очередь событий (class EventQueue) — основной структурный элемент, хранит как собственно очередь событий с временными метками их наступления (в секундах), так и указатель на соответствующий массив прямоугольников, а также его координаты на плоскости. Таким образом, очередь отвечает как за сохранение событий (функция push), за применение событий, произошедших к заданному моменту времени (функция next), так и за отрисовку связанного с ней массива прямоугольников (функция paint).

Зонд (class Probe) — сортируемый объект, генерирующий события с помощью перегруженных операторов. Событиям сравнения и присваивания нужны индексы элементов, поэтому зонд кроме собственно числового значения хранит и свой индекс. Кроме того, так как очередь событий требует «глобальную» временную метку, было решено перенести ответственность за вычисление времени наступления события на специальный класс TimeMapper, который получает объект IEvent, вычисляет время и уже сам отправляет событие в очередь. Поэтому зонд хранит указатель на объект TimeMapper, который является общим для всех элементов сортируемого массива.

Сцена (class Scene) — главный класс, соединяющий все куски визуализатора воедино. Хранит набор очередей событий и соответствующих массивов прямоугольников. Функция makeNew сбрасывает очереди, создаёт и сортирует новые массивы, заполняя очереди событиями. Функция paint отрисовывает все столбцы ( «всю сцену» ). Функция next передаёт время функциям next всех очередей, возвращает true, если хотя бы одно событие было извлечено из очередей (что используется как условие обновления буфера кадра в обработчике frameMove).

Таймер (class Timer) написан для удобного управления скоростью визуализации. Хранит «реальное» время, возвращаемое std::clock(), и «виртуальное» время визуализации в «секундах», с учётом скорости течения времени (поле factor).

UPDATE. 2013.01.23. Более новый вариант расположен здесь (в составе AV Orchid).

Графы-4: максимальный поток

Введение

Рассматривается ориентированный граф («сеть»), в котором есть две выделенные вершины: исток s и сток t. Исток не имеет входящих рёбер, сток — исходящих рёбер. Кроме того, если в сети существует ребро (u, v), то не существует ребра (v, u). Для каждого ребра (u, v) задана пропускная способность («capacity»), обозначаемая обычно c(u, v) или A(u, v). Ниже она задаётся матрицей, так же как и сам граф: пропускная способность несуществующих рёбер равна нулю.

Требуется найти максимальный поток F, который может пропустить сеть в силу заданных пропускных способностей. Поток в сети задаётся матрицей потоков, проходящих по рёбрам. Должно выполняться условие F(u, v) ≤ A(u, v).

Алгоритмы

Известно два основных класса алгоритмов, решающих задачу максимизации потока на сети:

  1. Класс алгоритмов Форда-Фалкерсона (Ford-Fulkerson algorithm), строящих дополняющие цепи.
  2. Алгоритмы проталкивания предпотока (preflow pushing algorithm, push-relabel algorithm), пытающихся максимально заполнить сеть, оперируя избытком «жидкости» в вершинах.

Алгоритмы Форда-Фалкерсона

Алгоритмы основаны на теореме Форда-Фалкерсона. Общий принцип работы:

  1. Устанавливается F(u, v) = 0, для всех пар u, v.
  2. Пока можно построить новую дополняющую цепь на остаточном графе, дополняем F(u, v) в силу построенной цепи.

Остаточный граф (residual graph) — граф с пропускными способностями рёбер равными A(u, v) — F(u, v), если A(u, v) > 0 (прямая дуга), или F(v, u), если A(u, v) = 0 (обратная дуга).

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

Разные реализации общей идеи по-разному решают вопрос поиска очередной дополняющей цепи. Существует два классических алгоритма 1970-х: алгоритм Эдмондса-Карпа (реализован в BGL) сложности O(N M2) ≤ O(N5) и алгоритм Ефима Диница сложности O(N2 M) ≤ O(N4). К новым алгоритмам этого семейства можно отнести алгоритм Бойкова-Колмогорова (описание из BGL, статья).

Алгоритмы проталкивания предпотока

В 1980-х появилось новое семейство алгоритмов построения максимального потока типа, отличного от ранее известных алгоритмов Форда-Фалкерсона. В данном контексте будем называть F предпотоком. Для предпотока выполняются следующие условия: F(u, v) ≤ A(u, v) и F(u, v) = —F(v, u).

Пусть для каждой вершины v задана высота h[v]  ≥ 0 и избыток e[v] ≥ 0. Вершины, в которых избыток больше нуля, называются активными. Для любых вершин u, v, кроме стока и истока, должно быть выполнено условие: если (u, v) принадлежит остаточному графу, то h[u] ≤ h[v] + 1. Обычно полагают h[s] = N, h[t] = 0. Избыток в вершине равен разностью между суммой потоков на исходящих рёбрах и суммой потоков на входящих рёбрах.

Основная идея алгоритма состоит в первоначальном максимальном заполнении исходящих рёбер истока и «проталкивании» избытка из активных вершин по рёбрам остаточного графа в вершины, имеющими высоту, меньшую на 1. Если весь избыток вывести не удалось, то высота вершины увеличивается. Алгоритм останавливается, когда не остаётся активных вершин.

Базовая версия алгоритма (реализация на e-maxx.ru) без упорядочивания активных вершин имеет оценку сложности O(N2 M) ≤ O(N4). Вариант, «прогоняющий» вершины через очередь имеет оценку сложности O(N3). Вариант, поднимающий в начало очереди активных вершин вершину, в случае увеличения её высоты (реализация на e-maxx.ru, описание реализации BGL), имеет оценку O(N2 sqrt(M)) ≤ O(N3).

Реализация

Была написана реализация алгоритмов Эдмондса-Карпа и проталкивания предпотока с переупорядочиванием поднятых вершин (maxFlowPushRelabel) и с обработкой активных вершин в порядке очереди («базовый», maxFlowPushing).

Алгоритмы реализованы как функции-шаблоны и имеют параметр шаблона WeightMatrix, задающий тип матрицы пропускных способностей.

Кроме того, алгоритм Форда-Фалкерсона принимает тип-параметр Augmentor, который считается классом, определяющим статическую функцию augment, строящую следующую дополняющую цепь. Имеется одна реализация: EdmondsKarp, строящая цепь поиском в ширину. Таким образом, функция maxFlowFordFulkerson<EdmondsKarp> реализует алгоритм Эдмондса-Карпа.

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

Задания

  1. Прочитать и понять код, почему первое выводимое число всегда 1, а последнее — 0?
  2. Каким ещё образом можно реализовать выборку наивысшей активной вершины?
  3. Можно ли объединить код maxFlowPushing и maxFlowPushRelabel в одну функцию-шаблон?

BY-1

скриншот

Введение

Результатом BY-0 стала простенькая программа-симулятор движения твёрдых сфер на плоскости. Очевидная следующая задача — сделать движение в трёхмерном пространстве. Сделаем это в два шага (1а и 1б).

Шаг 1а

Цель — реализовать движение в 3D. Будем сравнивать старый и новый варианты.

Первым делом расширим двумерный вектор Vec2 до трёхмерного вектора Vec3, добавив координату z. Ключевое слово explicit в определении конструктора предотвращает неявное преобразование скаляра к вектору. Теперь вызов

glTranslated(p.x, p.y, p.z);

вместо

glTranslated(p.x, p.y, 0.0);

в Ball::paint() позиционирует шар в трёхмерном пространстве.

Вместо действия одной центральной силы притяжения будем рассчитывать взаимное притяжение шаров друг к другу по формуле закона всемирного тяготения (естественно, значение константы G у нас своё). Сила, действующая на пару шаров вычисляется функцией gravity(Ball&, Ball&).

Кроме того, на шары будет действовать сила «вязкого» трения (не зависящая, впрочем, от взаимно-относительного положения тел) равная -Fr^2 v, где F — коэффициент трения, r — радиус шара, v — скорость шара. Сила трения вычисляется функцией friction(Ball&). Размеры окна в процессе вычисления сил нам более не нужны, поэтому переменные width и height убраны.

Рисование в 3D требует применения техник отсечения невидимых поверхностей. Стандартной техникой, работающей на попиксельном уровне и реализованной аппаратно в GPU, является тест глубины с применением буфера глубины (также называемого z-буфером). В OpenGL тест глубины включается с помощью команды (см. main())

glEnable(GL_DEPTH_TEST);

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

glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

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

  glEnable(GL_NORMALIZE);
  glEnable(GL_LIGHTING);
  glEnable(GL_LIGHT0);

Параметры источника света LIGHT0 задаются по умолчанию. Это направленный (бесконечно удалённый, заданный направлением лучей) источник белого цвета. Для изменения его параметров можно использовать функции glLight*. При расчёте освещения полагается, что нормали имеют единичную длину. Однако, например, при масштабировании объекта (мы задаём радиус шара), происходит масштабирование и нормалей в том числе, что изменяет их длину. Вызов glEnable(GL_NORMALIZE) предлагает OpenGL позаботиться об этом: включает автоматическую нормализацию нормалей.

Обработчик события изменения сторон окна теперь выставляет (пока всё ещё ортогональную) проекцию примерно фиксированного размера, отображая вдоль наибольшего измерения (ширины или высоты окна) координаты сцены от -500 до 500 (вдоль оси 0x или оси 0y, соответственно).

Что до прочих мелочей, то изменилась «политика выхода» из программы: вместо «жёсткого» exit(0) вызываем glutExit(). Кроме того, чтобы распределить вычислительную нагрузку на все доступные в системе аппаратно поддерживаемые потоки ЦП внутри функции frameMove() используется директива OpenMP:

# pragma omp parallel for

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

Шаг 1б

Следующая цель — управляемая проекция с имитацией перспективы, попросту говоря, «камера». Будем сравнивать вариант, полученный на шаге 1а, и окончательный вариант BY-1.

Камера (класс SphericalCoords) у нас будет довольно простая: основанная на сферических координатах. Положение её в пространстве задаётся двумя углами phi (долгота) и theta (широта), а также расстоянием до начала координат r. Пользователь может поворачивать камеру вокруг начала координат влево (A) и вправо (D) на угол dphi, вверх (Z) и вниз (X) на угол dtheta или придвигать (W) и отодвигать (S) её от начала координат домножением-делением r на заданный коэффициент rf.

Кроме клавиш, генерирующих при нажатии ASCII-код, можно использовать «специальные клавиши» (управление курсором, функциональные). Для этого требуется определить специальный обработчик (в нашей программе он назван specialKeyPressed) и подключить его к GLUT вызовом glutSpecialFunc. Теперь камерой можно управлять клавишами-стрелками.

Матрица проекции строится как произведение двух матриц.

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

gluPerspective(75, double(w) / double(h), 16.0, 4096.0);

который выставляет угол зрения по вертикали в 75 градусов, соотношение сторон (коэффициент, на который нужно домножить угол зрения по вертикали, чтобы получить угол зрения по горизонтали), расстояние до ближней отсекающей плоскости — 16 единиц (всё, что ближе к камере — не рисуется), расстояние до дальней отсекающей плоскости — 4096 единиц (всё, что дальше от камеры — не рисуется. Таким образом задаются размеры объёма отсечения, но не его положение и ориентация в пространстве. Изменяя угол зрения можно добиться эффекта оптического зума.

2. Матрицы камеры, задающей положение и ориентацию видимой области в пространстве. Эта матрица вычисляется в функции SphericalCoords::apply() при помощи функции gluLookAt, которая позволяет задать координаты положения камеры, координаты точки, на которую направлена (смотрит) камера и направление «вверх» (с помощью выбора которого можно вращать камеру вокруг оси зрения). Как нетрудно догадаться, фактически эта функция вычисляет матрицу поворота и сдвига.

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

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

glEnable(GL_BLEND);

Формула смешивания задаётся с помощью функции glBlendFunc:

glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);

Первый параметр определяет, на что будет домножен цвет «источника» (нового пикселя), в данном случае — на собственное значение альфа-канала. Второй параметр определяет, на что будет домножен цвет «назначения» (пикселя в буфере кадра), в данном случае — на (1 — a), где a — значение альфа-канала нового пикселя. Цвет результирующего пикселя получается как сумма этих двух произведений. Например, вызов

glBlendFunc(GL_ONE, GL_ONE);

приведёт к простому сложению цветов пикселей.

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

glEnable(GL_CULL_FACE);

Таким образом, полигоны, нормали к которым «смотрят» от зрителя (эти полигоны повёрнуты задней стороной), не будут отрисовываться, что уберёт артефакты при рисовании одиночных шаров (однако не решит проблем с наложением нескольких шаров).

BY-0

Введение

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

Шарики

Изначально код рассчитан на N шариков. При этом почему-то было сделано так, что столкновения обрабатываются только между первыми двумя. Впрочем, это несложно исправить.

Вначале генерируется случайная «целевая» точка столкновения внутри окна (функция new_circles). N шариков со случайными массами и радиусами расставляются в случайных позициях со скоростями, направленными в точку столкновения, по величине равными четверти расстояния от центра шарика до точки столкновения. Таким образом, шарики должны столкнуться примерно через 4 секунды.

За рисование (рендеринг) отвечает функция repaint. Собственно же движения и столкновения шариков рассчитываются в функции frame_move. Время на момент расчёта последнего кадра хранится в глобальной переменной t. В случае, если на момент вызова frame_move() прошло больше 0 микросекунд, происходит изменение положения шариков соответственно их скоростям и прошедшему времени. Если центры шариков 0 и 1 находятся друг от друга не дальше суммы радиусов этих шариков, то регистрируется столкновение и скорости шариков меняются соответственно. Код под комментарием «improve distance» производит «разведение» шариков, чтобы они перестали касаться. Работа этого кода может приводить к «туннелированию» шариков, если результирующие скорости оказались близки к параллельным.

Функция handle_key реализует реакцию на нажатие пользователем клавиш, когда активно окно, где происходит рисование шариков (под MS Visual Studio данный файл следует компилировать как консольное приложение, соответственно будет два окна: консоль и окно с графикой, правда, чтобы его скомпилировать и успешно запустить, нужно раздобыть библиотеку GLUT в составе glut.h, glut32.lib и glut32.dll). Нажатие на enter или пробел создаёт новый набор шариков, нажатие на escape (код 27 в ASCII) или q приводит к выходу из программы.

Функция main производит инициализацию GLUT, состояния контекста OpenGL, создаёт аппроксимацию сферы (единичного радиуса, 12 вершин на «круг» по долготе, 4 таких круга по широте), регистрирует (передаёт библиотеке указатели на эти функции) функции-обработчики событий (callbacks; под комментарием «register callbacks» ). Вызов glutMainLoop() запускает цикл-обработчик событий, из которого, собственно, и вызываются зарегистрированные нами функции.

Окна

Зачем нам нужен GLUT? Что до OpenGL, то это низкоуровневый графический API, обращающийся непосредственно к драйверу GPU. В его «компетенции» не входит создание окон GUI или, тем более, организация взаимодействия с пользователем (реакция на события клавиатуры и мыши). Всё это, а также инициализация виртуального устройства OpenGL (создание контекста OpenGL) реализуется операционной системой. Если же хочется обеспечить независимость от ОС и/или простоту кода (что важно для маленьких приложений), имеет смысл воспользоваться сторонними библиотеками-обёртками, предоставляющими указанный функционал. Например, библиотекой GLUT.

Библиотека GLUT имеет давно устоявшийся (в настоящее время не развивается), достаточно простой в использовании (самый простой?) интерфейс на языке C (можно использовать в программах на разных языках программирования). Конечно, есть другие библиотеки, с помощью которых можно сделать то же самое и что-нибудь ещё, как маленькие и простые (например, SFML), так и огромные, предоставляющие широкий функционал для создания серьёзных GUI приложений (например, Qt).

Новый вариант

Итак, старый вариант был переправлен.

Для сравнения старого и нового вариантов удобно использовать какой-нибудь вариант утилиты diff (у Microsoft есть своя реализация — WinDiff, которая может присутствовать в составе Visual Studio/Windows SDK; мне же больше нравится WinMerge).

В качестве реализации GLUT было решено взять более современную (чем оригинальная) open-source библиотеку freeglut (тем более, что есть сборка под Visual C++). Вместо микросекундного времени boost::date_time, будем использовать обычный миллисекундный std::clock.

Макрос GLUTCALLBACK нам теперь без надобности (был на случай особого соглашения вызова для функций-обработчиков, например, __stdcall).

Посмотрим Ball::paint vs circle::paint. Стек матриц не используется, поэтому перед установкой матрицы для шара, просто заменяем текущую матрицу на единичную. Далее, было:

glScaled(r, r, r);
glTranslated(p.x / r, p.y / r, 0.0);

стало

glTranslated(p.x, p.y, 0.0);
glScaled(r, r, r);

Указанные функции домножают текущую матрицу на матрицу соответствующего преобразования слева. Поэтому (новый вариант), сначала промасштабируем единичную сферу, сделав из неё сферу радиуса r, потом переместим её из начала координат в позицию (p.x, p.y), а не наоборот, как было в старом варианте (поэтому там деление на r).

Появилась новая функция Ball::move, которая осуществляет расчёт сдвига шара в соответствии с заданным ускорением a. Ускорение считается постоянным в течении указанного промежутка времени time_step. Метод первого порядка, его точность низкая, но вполне соответствует обычной «игровой физике». Итак, сначала нужно записать в a сумму ускорений, порождаемых силами, действующими на шар, потом вызвать move.

Сила здесь будет одна — похожее на гравитацию (прямо пропорциональное массе шара и константе G, обратно пропорциональное квадрату расстояния) притяжение к центру экрана. Вычислением ускорения, порождаемого этой силой, занимается функция applyCentralForce.

Функция frame_move переименована в frameMove и теперь рассчитывает столкновения для всех пар шариков. Кроме того, вместо того, чтобы «разводить» шарики, «зашедшие» внутрь друг друга, фиксация столкновения производится только в случае сближения центров шариков. Шарики могут оказаться внутри друг друга (что особенно заметно из-за того, что при случайном выборе позиции факт «попадания» внутрь не проверяется), зато нет «туннелирования».

Добавлен обработчик события изменения размеров окна рендеринга reshape, именно там устанавливается матрица проекции, которая в старом варианте просто не использовалась. Текущие размеры окна хранятся в глобальных переменных width и height.

С фактом наличия глобальных переменных пока решено смириться.