Реализация алгоритма двумерной упаковки Skyline

Упаковка 2D-прямоугольников в прямоугольники большего фиксированного размера необходима в большинстве мультимедийных проектов. В программировании GPU изменение текстур (binding) — затратный процесс. Поэтому при рендеринге текста не стоит использовать по одной текстуре на глиф, вместо этого желательно упаковать глифы в единую текстуру, называемую атласом. В 2D-играх содержащие спрайты атласы называются листами спрайтов (spritesheet). Листы спрайтов также используются для веб-сайтов, потому что скачивать один большой файл удобнее, чем по одному файлу на каждый значок/логотип.

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

Самым ценным ресурсом из найденных мной стал превосходный обзор Юкки Йулянки. В нём описано четыре типа алгоритмов и их практическая оценка. Выделяются из них два:

• MAXRECTS если вы знаете заранее, какие прямоугольники будете упаковывать («офлайн-упаковка»)

• SKYLINE если не знаете («онлайн-упаковка»)

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

Именно поэтому я остановился на алгоритме skyline. Он применяется в stb_rect_pack.h, fontstash, а значит, и в nanovg.

В этой статье объяснён алгоритм skyline и представлена его реализация. Реализация доступна онлайн и в общественном достоянии (UNLICENSE).

API, который нужно реализовать

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

struct JvPack2d {     // Задаём массив, способный хранить '2 * maxWidth' uint16_t     uint16_t* pSkyline;      // Размер атласа.     uint16_t maxWidth;     uint16_t maxHeight;      // Приватное состояние, инициализируемое нулём.     bool _bInitialized;     uint16_t _skylineCount; };  bool jvPack2dAdd(JvPack2d* pP2D, uint16_t const size[2], uint16_t pos[2]);

С maxWidth и maxHeight всё просто: это размеры целевого прямоугольника, например текстуры. Я решил использовать координаты в uint16_t, потому что хочу работать с GPU-текстурами, размер которых ограничен 16384 пикселями по каждой из сторон.

Массив pSkyline возлагает на пользователя ответственность за распределение. Требуемый объём 2 * maxWidth — это наихудший сценарий, при котором все прямоугольники имеют длину и ширину 1 пиксель. Если вы хотите сэкономить немного килобайтов и знаете, что все ширины прямоугольников будут кратны четырём (или можете это обеспечить), то можно разделить все ширины и maxWidth на 4 и умножить на 4 горизонтальные позиции, задаваемые jvPack2dAdd().

Как работает алгоритм Skyline

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

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

Алгоритм можно разделить на два этапа:

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

2. Обновляем структуру данных skyline силуэтом нового прямоугольника.

Структура данных

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

На языке C этот массив представлен буфером pSkyline, содержащим попеременно позиции X и Y точек, и количеством хранимых в буфере точек: skylineCount. По горизонтальной оси ширина maxWidth пикселей, то есть линия горизонта может иметь максимум maxWidth отрезков, если каждый отрезок имеет длину 1 пиксель, например, в форме зубчатых стен замка. Таким образом, массив линии горизонта содержит не более maxWidth двухмерных координат. Именно поэтому в API pSkyline должен иметь возможность хранить 2 * maxWidth чисел.

bool jvPack2dAdd(JvPack2d* pP2D, uint16_t const size[2], uint16_t pos[2]) {     // Не имеет смысла упаковывать прямоугольники с шириной/высотой 0, выполняем ранний выход.     uint16_t width = size[0];     uint16_t height = size[1];     if (width == 0 || height == 0)         return false;      // Вспомогательная структура, повышающая читаемость индексирования массива.     typedef struct Point {         uint16_t x;         uint16_t y;     } Point;     Point* pSkyline = (Point*)pP2D->pSkyline;      // Изначальное состояние содержит нижнюю левую точку.     if (!pP2D->_bInitialized) {         _jvASSERT(pP2D->maxWidth, >, 0);         _jvASSERT(pP2D->maxHeight, >, 0);         pP2D->_bInitialized = true;         pP2D->_skylineCount = 1;         pSkyline[0].x = 0;         pSkyline[0].y = 0;     }

Поиск точки вставки

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

Взяв прямоугольник размером (w,h), который пытаемся уместить в каждую возможную позицию (x,y), мы находим точки линии горизонта (x₂​,y_2​), которые горизонтально пересекаются с прямоугольником, то есть x≤x₂​<x+w, по двум причинам:

1. Если y₂​>y, то у прямоугольника есть коллизия с линией горизонта. В этом случае прямоугольник поднимается выше, чтобы предотвратить коллизию y←y₂​.

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

В коде интервал накладывающихся точек для лучшего на текущий момент кандидата хранится как два индекса idxBest (включая его) и idxBest2 (исключая его). Эти накладывающиеся точки всегда представляют последовательный интервал, потому что массив отсортирован по увеличению горизонтальной позиции.

    // Псевдонимы делают код более кратким.     uint16_t maxWidth = pP2D->maxWidth;     uint16_t maxHeight = pP2D->maxHeight;     uint16_t skylineCount = pP2D->_skylineCount;      // Хранит наилучшего на текущий момент кандидата.     uint16_t idxBest = UINT16_MAX;     uint16_t idxBest2 = UINT16_MAX;     uint16_t bestX = UINT16_MAX;     uint16_t bestY = UINT16_MAX;      // Цикл поиска наилучшего места.     for (uint16_t idx = 0; idx < skylineCount; ++idx) {         uint16_t x = pSkyline[idx].x;         uint16_t y = pSkyline[idx].y;          if (width > maxWidth - x)             break; // Мы достигли правой границы атласа.         if (y >= bestY)             continue; // Нам не найти кандидата лучше, чем текущий.          // Поднимаем 'y' так, чтобы находиться выше всех точек пересечений.         uint16_t xMax = x + width;         uint16_t idx2;         for (idx2 = idx + 1; idx2 < skylineCount; ++idx2) {             if (xMax <= pSkyline[idx2].x)                 break; // Мы не достигли следующих точек.             if (y < pSkyline[idx2].y)                 y = pSkyline[idx2].y; // Поднимаем 'y' так, чтобы не было пересечений.         }          if (y >= bestY)             continue; // Мы не нашли кандидата лучше, чем текущий.         if (height > maxHeight - y)             continue; // Мы достигли верхней границы атласа.          idxBest = idx;         idxBest2 = idx2;         bestX = x;         bestY = y;     }      if (idxBest == UINT16_MAX)         return false; // Не нашли свободного места.     _jvASSERT(idxBest, <, idxBest2);     _jvASSERT(idxBest2, >, 0);

Обновление структуры данных линии горизонта

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

1. Новый прямоугольник C накладывается на две точки, которые будут удалены. Добавляются две точки: TL (левая верхняя) и BR (правая нижняя). Обратите внимание, что BR опущена на тот же ряд, что и последняя накладывающаяся точка.

2. Новый прямоугольник D накладывается на одну позицию (точку начала прямоугольника), которая будет удалена. Добавляется левая верхняя точка TL. Правая нижняя точка не добавляется BR, потому что она будет соответствовать правой границе атласа.

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

Выбор решения соответствует следующему коду.

    // Мы заменяем точки, затенённые текущим rect, на новые точки.      uint16_t removedCount = idxBest2 - idxBest;      Point newTL, newBR; // TopLeft, BottomRight     newTL.x = bestX;     newTL.y = bestY + height;     newBR.x = bestX + width;     newBR.y = pSkyline[idxBest2 - 1].y;     bool bBottomRightPoint =         (idxBest2 < skylineCount ? newBR.x < pSkyline[idxBest2].x : newBR.x < maxWidth);     // TopLeft вставляется всегда     uint16_t insertedCount = 1 + bBottomRightPoint;      _jvASSERT(skylineCount + insertedCount - removedCount, <=, maxWidth);

Затем нам нужно или уменьшить, или расширить массив, чтобы уместить удаление и вставку. Так как в стандартной библиотеке есть средства манипуляций с массивами наподобие std::vector из C++, это может быть одним простым вызовом .erase() с последующим .insert(). Но мне хочется, чтобы в этой статье алгоритм был выражен вне зависимости от стандартной библиотеки, поэтому реализация сдвигает элементы явным образом:

    // Логика вставки и удаления в массиве.     // Если бы мы использовали стандартную библиотеку, это можно было бы сделать в двух строках,     // например, в C++ это были бы std::vector erase() + insert(),     // но я предпочитаю, чтобы алгоритм был явным и автономным.      if (insertedCount > removedCount) {         // Расширение. Сдвигаем элементы вправо. Нам нужно выполнять итерации в обратную сторону.         uint16_t idx = skylineCount - 1;         uint16_t idx2 = idx + (insertedCount - removedCount);         for (; idx >= idxBest2; --idx, --idx2)             pSkyline[idx2] = pSkyline[idx];         pP2D->_skylineCount = skylineCount + (insertedCount - removedCount);     }     else if (insertedCount < removedCount) {         // Уменьшение. Сдвигаем элементы влево. нам нужно выполнять итерации вперёд.         uint16_t idx = idxBest2;         uint16_t idx2 = idx - (removedCount - insertedCount);         for (; idx < skylineCount; ++idx, ++idx2)             pSkyline[idx2] = pSkyline[idx];         pP2D->_skylineCount = skylineCount - (removedCount - insertedCount);     }      pSkyline[idxBest] = newTL;     if (bBottomRightPoint)         pSkyline[idxBest + 1] = newBR;      pos[0] = bestX;     pos[1] = bestY;     return true; }

Результаты

Ниже показана GIF-анимация решения алгоритмом skyline исходной задачи:

У меня нет подходящего бенчмарка для измерения производительности, но в рамках юнит-тестов я написал четыре «наихудших» сценария, в которых мы заполняем атлас различными паттернами, и из них было достаточно легко получить измерения. Строго говоря, там измеряются и утверждения (assertion), но после профилирования кода выяснилось, что assertion не находятся на «горячих» путях исполнения кода, поэтому измерения как будто похожи на правду. Сценарии параметризируются одним параметром ATLAS_SIZE.

1. WorstCaseWidth: размер атласа равен (ATLAS_SIZE, 2), и мы вставляем прямоугольники размера (1,1). Упаковывается 2 * ATLAS_SIZE прямоугольников. Требуется 4 * ATLAS_SIZE байтов.

2. WorstCaseHeight: размер атласа равен (2, ATLAS_SIZE), и мы вставляем прямоугольники размера (1,1). Упаковывается 2 * ATLAS_SIZE прямоугольников. Требуется 8 байтов.

3. WorstCaseDiagonalVertical: размер атласа равен (ATLAS_SIZE, ATLAS_SIZE), и мы вставляем прямоугольники шириной 1 пиксель для создания диагонального паттерна. Упаковывается 2 * ATLAS_SIZE - 1 прямоугольников. Требуется 4 * ATLAS_SIZE байтов.

4. WorstCaseDiagonalHorizontal: размер атласа — (ATLAS_SIZE, ATLAS_SIZE), и мы вставляем прямоугольники шириной 1 пиксель для создания диагонального паттерна. Упаковывается 2 * ATLAS_SIZE - 1 прямоугольников. Требуется 4 * ATLAS_SIZE байтов.

Вот сравнение времени выполнения соответствующих юнит-тестов. С каждой строкой мы удваиваем ATLAS_SIZE. То есть количество пакуемых прямоугольников тоже удваивается. Измерения выполнялись на ноутбуке с процессором AMD Ryzen 5 7520U 2.80 GHz, то есть сырое время не так релевантно, но эволюция с увеличением ATLAS_SIZE любопытна.

• WorstCaseHeight и WorstCaseDiagonalHorizontal составлены из двух прямоугольников на строку. Вставки выполняются быстро и время линейно зависит от ATLAS_SIZE; затраты на каждый добавляемый прямоугольник фиксированы.

• WorstCaseWidth и WorstCaseDiagonalVertical составлены из двух прямоугольников на столбец. Время вставки гораздо выше и увеличивается квадратично от ATLAS_SIZE; затраты на каждый добавляемый прямоугольник линейны относительно ATLAS_SIZE.

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

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